{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Quantitative Finance\n",
    "\n",
    "Copyright (c) 2019 Python Charmers Pty Ltd, Australia, <https://pythoncharmers.com>. All rights reserved.\n",
    "\n",
    "<img src=\"img/python_charmers_logo.png\" width=\"300\" alt=\"Python Charmers Logo\">\n",
    "\n",
    "Published under the Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) license. See `LICENSE.md` for details.\n",
    "\n",
    "Sponsored by Tibra Global Services, <https://tibra.com>\n",
    "\n",
    "<img src=\"img/tibra_logo.png\" width=\"300\" alt=\"Tibra Logo\">\n",
    "\n",
    "\n",
    "## Module 1.3: Ordinary Least Squares\n",
    "\n",
    "### 1.3.1 Ordinary Least Squares\n",
    "\n",
    "The Ordinary Least Squares (hereafter OLS) algorithm is a key algorithm for prediction, regression and many other learning tasks. Despite it initially looking like *just* a linear algorithm, it can be extended to be able to predict a wide variety of complex tasks. It is a powerful algorithm with some solid mathematical proofs behind it, and is a very practically useful algorithm for prediction and modelling tasks.\n",
    "\n",
    "We will review the OLS algorithm in significant depth over the next few notebooks, and use it extensively over the rest of this course.\n",
    "The OLS algorithm aims to solve problems of the form:\n",
    "\n",
    "$ Y = X\\beta + \\boldsymbol{u}$\n",
    "\n",
    "Where:\n",
    "* $Y$ is the variable we are attempting to predict (the dependent variable)\n",
    "* $X$ is our independent variable (or multiple variables, as we will see in the next notebook)\n",
    "* $\\beta$ is the parameters we are trying to learn\n",
    "* $u$ is the error of our model, as we can often not predict data perfectly.\n",
    "\n",
    "This type of model is a **Linear Regression Model**, as the model is a linear combination of $X$ and $\\beta$. Note, however, that the values of $X$ do not have to be linear, for instance, we can used $x^2$ as our $X$ value, and we still use Linear Regression to learn that model (although in some cases it may no longer be the best algorithm - a problem for later). We use OLS to find the best values for $\\beta$ for our Linear Regression Model.\n",
    "\n",
    "OLS aims to learn the *best* values for $\\beta$, which is defined as the values that minimise the value $u$, our error. Before we get to how to do that, we will have a look at the Linear Regression Model in more detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression Model\n",
    "\n",
    "In a Linear Regression Model, we have our independent $X$ values, and our dependent $Y$ values, and we wish to find a linear relationship between them. For instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run setup.ipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "prices = pd.read_hdf(\"data/HousePricesVsInterestRates.hdf\", key=\"NSW\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>MedianPrice</th>\n",
       "      <th>CashRate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2007-03-31</th>\n",
       "      <td>517447.312500</td>\n",
       "      <td>6.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2007-06-30</th>\n",
       "      <td>528379.733333</td>\n",
       "      <td>6.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2007-09-30</th>\n",
       "      <td>511257.021277</td>\n",
       "      <td>6.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2007-12-31</th>\n",
       "      <td>547164.071090</td>\n",
       "      <td>6.75</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2008-03-31</th>\n",
       "      <td>567444.245132</td>\n",
       "      <td>7.25</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              MedianPrice  CashRate\n",
       "2007-03-31  517447.312500      6.25\n",
       "2007-06-30  528379.733333      6.25\n",
       "2007-09-30  511257.021277      6.50\n",
       "2007-12-31  547164.071090      6.75\n",
       "2008-03-31  567444.245132      7.25"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prices.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Divide the MedianPrice by 1000 to get easier numbers to work with\n",
    "prices['MedianPrice'] /= 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "  #altair-viz-01a90dcee8a642d48bd5cc9e1eb5121c.vega-embed {\n",
       "    width: 100%;\n",
       "    display: flex;\n",
       "  }\n",
       "\n",
       "  #altair-viz-01a90dcee8a642d48bd5cc9e1eb5121c.vega-embed details,\n",
       "  #altair-viz-01a90dcee8a642d48bd5cc9e1eb5121c.vega-embed details summary {\n",
       "    position: relative;\n",
       "  }\n",
       "</style>\n",
       "<div id=\"altair-viz-01a90dcee8a642d48bd5cc9e1eb5121c\"></div>\n",
       "<script type=\"text/javascript\">\n",
       "  var VEGA_DEBUG = (typeof VEGA_DEBUG == \"undefined\") ? {} : VEGA_DEBUG;\n",
       "  (function(spec, embedOpt){\n",
       "    let outputDiv = document.currentScript.previousElementSibling;\n",
       "    if (outputDiv.id !== \"altair-viz-01a90dcee8a642d48bd5cc9e1eb5121c\") {\n",
       "      outputDiv = document.getElementById(\"altair-viz-01a90dcee8a642d48bd5cc9e1eb5121c\");\n",
       "    }\n",
       "    const paths = {\n",
       "      \"vega\": \"https://cdn.jsdelivr.net/npm/vega@5?noext\",\n",
       "      \"vega-lib\": \"https://cdn.jsdelivr.net/npm/vega-lib?noext\",\n",
       "      \"vega-lite\": \"https://cdn.jsdelivr.net/npm/vega-lite@5.16.3?noext\",\n",
       "      \"vega-embed\": \"https://cdn.jsdelivr.net/npm/vega-embed@6?noext\",\n",
       "    };\n",
       "\n",
       "    function maybeLoadScript(lib, version) {\n",
       "      var key = `${lib.replace(\"-\", \"\")}_version`;\n",
       "      return (VEGA_DEBUG[key] == version) ?\n",
       "        Promise.resolve(paths[lib]) :\n",
       "        new Promise(function(resolve, reject) {\n",
       "          var s = document.createElement('script');\n",
       "          document.getElementsByTagName(\"head\")[0].appendChild(s);\n",
       "          s.async = true;\n",
       "          s.onload = () => {\n",
       "            VEGA_DEBUG[key] = version;\n",
       "            return resolve(paths[lib]);\n",
       "          };\n",
       "          s.onerror = () => reject(`Error loading script: ${paths[lib]}`);\n",
       "          s.src = paths[lib];\n",
       "        });\n",
       "    }\n",
       "\n",
       "    function showError(err) {\n",
       "      outputDiv.innerHTML = `<div class=\"error\" style=\"color:red;\">${err}</div>`;\n",
       "      throw err;\n",
       "    }\n",
       "\n",
       "    function displayChart(vegaEmbed) {\n",
       "      vegaEmbed(outputDiv, spec, embedOpt)\n",
       "        .catch(err => showError(`Javascript Error: ${err.message}<br>This usually means there's a typo in your chart specification. See the javascript console for the full traceback.`));\n",
       "    }\n",
       "\n",
       "    if(typeof define === \"function\" && define.amd) {\n",
       "      requirejs.config({paths});\n",
       "      require([\"vega-embed\"], displayChart, err => showError(`Error loading script: ${err.message}`));\n",
       "    } else {\n",
       "      maybeLoadScript(\"vega\", \"5\")\n",
       "        .then(() => maybeLoadScript(\"vega-lite\", \"5.16.3\"))\n",
       "        .then(() => maybeLoadScript(\"vega-embed\", \"6\"))\n",
       "        .catch(showError)\n",
       "        .then(() => displayChart(vegaEmbed));\n",
       "    }\n",
       "  })({\"config\": {\"view\": {\"continuousWidth\": 300, \"continuousHeight\": 300}}, \"data\": {\"url\": \"altair-data-16a52107771e16dd783182898e95326f.json\", \"format\": {\"type\": \"json\"}}, \"mark\": {\"type\": \"circle\"}, \"encoding\": {\"x\": {\"field\": \"CashRate\", \"type\": \"quantitative\"}, \"y\": {\"field\": \"MedianPrice\", \"type\": \"quantitative\"}}, \"height\": 400, \"width\": 800, \"$schema\": \"https://vega.github.io/schema/vega-lite/v5.16.3.json\"}, {\"mode\": \"vega-lite\"});\n",
       "</script>"
      ],
      "text/plain": [
       "alt.Chart(...)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alt.Chart(prices).mark_circle().encode(\n",
    "    x='CashRate',\n",
    "    y='MedianPrice'\n",
    ").properties(\n",
    "    width=800,\n",
    "    height=400\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this data, there is a linear relationship between the Cash Interest Rate and House Prices that we will attempt to model.\n",
    "\n",
    "Side note: while there is a causual relationship between the two (a lower cash rate gives more people the ability to take out loans), the linear relationship here is a drastic oversimplification of things. For instance, the interest rate has been in relatively steady decline for a while, while house prices have gone up, and generally would have due to inflation anyway. A more formal analysis would need to remove the effect of inflation and other impacts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "  #altair-viz-a0bd35cd50de4af09cbf7ac5c0e5df1c.vega-embed {\n",
       "    width: 100%;\n",
       "    display: flex;\n",
       "  }\n",
       "\n",
       "  #altair-viz-a0bd35cd50de4af09cbf7ac5c0e5df1c.vega-embed details,\n",
       "  #altair-viz-a0bd35cd50de4af09cbf7ac5c0e5df1c.vega-embed details summary {\n",
       "    position: relative;\n",
       "  }\n",
       "</style>\n",
       "<div id=\"altair-viz-a0bd35cd50de4af09cbf7ac5c0e5df1c\"></div>\n",
       "<script type=\"text/javascript\">\n",
       "  var VEGA_DEBUG = (typeof VEGA_DEBUG == \"undefined\") ? {} : VEGA_DEBUG;\n",
       "  (function(spec, embedOpt){\n",
       "    let outputDiv = document.currentScript.previousElementSibling;\n",
       "    if (outputDiv.id !== \"altair-viz-a0bd35cd50de4af09cbf7ac5c0e5df1c\") {\n",
       "      outputDiv = document.getElementById(\"altair-viz-a0bd35cd50de4af09cbf7ac5c0e5df1c\");\n",
       "    }\n",
       "    const paths = {\n",
       "      \"vega\": \"https://cdn.jsdelivr.net/npm/vega@5?noext\",\n",
       "      \"vega-lib\": \"https://cdn.jsdelivr.net/npm/vega-lib?noext\",\n",
       "      \"vega-lite\": \"https://cdn.jsdelivr.net/npm/vega-lite@5.16.3?noext\",\n",
       "      \"vega-embed\": \"https://cdn.jsdelivr.net/npm/vega-embed@6?noext\",\n",
       "    };\n",
       "\n",
       "    function maybeLoadScript(lib, version) {\n",
       "      var key = `${lib.replace(\"-\", \"\")}_version`;\n",
       "      return (VEGA_DEBUG[key] == version) ?\n",
       "        Promise.resolve(paths[lib]) :\n",
       "        new Promise(function(resolve, reject) {\n",
       "          var s = document.createElement('script');\n",
       "          document.getElementsByTagName(\"head\")[0].appendChild(s);\n",
       "          s.async = true;\n",
       "          s.onload = () => {\n",
       "            VEGA_DEBUG[key] = version;\n",
       "            return resolve(paths[lib]);\n",
       "          };\n",
       "          s.onerror = () => reject(`Error loading script: ${paths[lib]}`);\n",
       "          s.src = paths[lib];\n",
       "        });\n",
       "    }\n",
       "\n",
       "    function showError(err) {\n",
       "      outputDiv.innerHTML = `<div class=\"error\" style=\"color:red;\">${err}</div>`;\n",
       "      throw err;\n",
       "    }\n",
       "\n",
       "    function displayChart(vegaEmbed) {\n",
       "      vegaEmbed(outputDiv, spec, embedOpt)\n",
       "        .catch(err => showError(`Javascript Error: ${err.message}<br>This usually means there's a typo in your chart specification. See the javascript console for the full traceback.`));\n",
       "    }\n",
       "\n",
       "    if(typeof define === \"function\" && define.amd) {\n",
       "      requirejs.config({paths});\n",
       "      require([\"vega-embed\"], displayChart, err => showError(`Error loading script: ${err.message}`));\n",
       "    } else {\n",
       "      maybeLoadScript(\"vega\", \"5\")\n",
       "        .then(() => maybeLoadScript(\"vega-lite\", \"5.16.3\"))\n",
       "        .then(() => maybeLoadScript(\"vega-embed\", \"6\"))\n",
       "        .catch(showError)\n",
       "        .then(() => displayChart(vegaEmbed));\n",
       "    }\n",
       "  })({\"config\": {\"view\": {\"continuousWidth\": 300, \"continuousHeight\": 300}}, \"layer\": [{\"data\": {\"url\": \"altair-data-16a52107771e16dd783182898e95326f.json\", \"format\": {\"type\": \"json\"}}, \"mark\": {\"type\": \"circle\"}, \"encoding\": {\"x\": {\"field\": \"CashRate\", \"type\": \"quantitative\"}, \"y\": {\"field\": \"MedianPrice\", \"type\": \"quantitative\"}}}, {\"data\": {\"url\": \"altair-data-0983ca0ffb151542a56784fdf2d68f3b.json\", \"format\": {\"type\": \"json\"}}, \"mark\": {\"type\": \"line\"}, \"encoding\": {\"x\": {\"field\": \"x\", \"type\": \"quantitative\"}, \"y\": {\"field\": \"y\", \"type\": \"quantitative\"}}}], \"height\": 300, \"width\": 800, \"$schema\": \"https://vega.github.io/schema/vega-lite/v5.16.3.json\"}, {\"mode\": \"vega-lite\"});\n",
       "</script>"
      ],
      "text/plain": [
       "alt.LayerChart(...)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_chart = alt.Chart(prices).mark_circle().encode(\n",
    "    x='CashRate',\n",
    "    y='MedianPrice'\n",
    ")\n",
    "\n",
    "\n",
    "x = np.arange(0, 8.25, 0.25)\n",
    "\n",
    "m = -100\n",
    "c = 1100\n",
    "\n",
    "y = m * x + c\n",
    "\n",
    "model = pd.DataFrame({'x': x, 'y':y})\n",
    "\n",
    "\n",
    "linear_model_chart = alt.Chart(model).mark_line().encode(\n",
    "    x='x', y='y'\n",
    ")\n",
    "\n",
    "(data_chart + linear_model_chart).properties(\n",
    "    width=800,\n",
    "    height=300\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quite a good fit! We can then work out our error, by using our \"model\" where $\\beta=4$ and a constant of 300 (we will see how this fits into a model soon), on the actual hours studied, and compare that to the actual scores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The error is 4935.24\n"
     ]
    }
   ],
   "source": [
    "predicted_prices = prices['CashRate'] * m + c  # Same values as before - this is our \"model\"\n",
    "error = np.sum(np.abs(predicted_prices - prices['MedianPrice']))  # The \"absolute error\"\n",
    "print(\"The error is {:.2f}\".format(error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The OLS algorithm aims to minimise the squared error, which is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The squared error is 689865.92\n"
     ]
    }
   ],
   "source": [
    "squared_error = np.sum((predicted_prices - prices['MedianPrice']) ** 2)\n",
    "print(\"The squared error is {:.2f}\".format(squared_error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both errors are valid to use, but the squared error is nicer from an algebraic perspective, allowing for gradients to be computed, and automatically account for positive/negative errors in the same way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercises\n",
    "\n",
    "1. Why wouldn't you use the actual error, that is `np.sum(predicted_prices - prices['MedianPrice'])`?\n",
    "2. Try three different models (change the values of m and c in the equation before). What is your best performing model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The error is 4935.24\n",
      "The squared error is 689865.92\n",
      "The error is 5181.09\n",
      "The squared error is 906373.12\n",
      "The error is 4083.08\n",
      "The squared error is 468521.42\n"
     ]
    }
   ],
   "source": [
    "m = -100\n",
    "c = 1100\n",
    "predicted_prices = prices['CashRate'] * m + c  # Same values as before - this is our \"model\"\n",
    "error = np.sum(np.abs(predicted_prices - prices['MedianPrice']))  # The \"absolute error\"\n",
    "print(\"The error is {:.2f}\".format(error))\n",
    "\n",
    "squared_error = np.sum((predicted_prices - prices['MedianPrice']) ** 2)\n",
    "print(\"The squared error is {:.2f}\".format(squared_error))\n",
    "\n",
    "m = -110\n",
    "c = 1000\n",
    "predicted_prices = prices['CashRate'] * m + c  # Same values as before - this is our \"model\"\n",
    "error = np.sum(np.abs(predicted_prices - prices['MedianPrice']))  # The \"absolute error\"\n",
    "print(\"The error is {:.2f}\".format(error))\n",
    "\n",
    "squared_error = np.sum((predicted_prices - prices['MedianPrice']) ** 2)\n",
    "print(\"The squared error is {:.2f}\".format(squared_error))\n",
    "\n",
    "m = -90\n",
    "c = 1000\n",
    "predicted_prices = prices['CashRate'] * m + c  # Same values as before - this is our \"model\"\n",
    "error = np.sum(np.abs(predicted_prices - prices['MedianPrice']))  # The \"absolute error\"\n",
    "print(\"The error is {:.2f}\".format(error))\n",
    "\n",
    "squared_error = np.sum((predicted_prices - prices['MedianPrice']) ** 2)\n",
    "print(\"The squared error is {:.2f}\".format(squared_error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/linear_regression_intro.py`*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ordinary Least Squares Derivation\n",
    "\n",
    "Here we will go through the derivation of the OLS algorithm with matrix mathematics. We will extend this further in the next notebook where we look at multivariate OLS. Before we do that, we will extend our Linear Regression Model slightly to account for our constant value (the \"c\" in our model for house prices):\n",
    "\n",
    "$ Y = X\\beta + \\boldsymbol{u}$  (as before)\n",
    "\n",
    "where $\\beta = [\\beta_1, \\beta_2]$ and $X = [1, X]$.\n",
    "\n",
    "In other words, $\\beta$ is two values - the multiplier against the constant (hereafter just \"the constant\") and the value that the independent variable is multiplied against ($X\\beta_2$). We extend $X$ by adding a column of 1s to it, and then use matrix multiplication, as we will review in the next notebook. This multiplies the 1s by the constant and also the independent variable by its coefficient.\n",
    "\n",
    "The aim of OLS is to choose $\\beta$ to minimise $u^2$ (the squared error), which we can do via the following equations:\n",
    "\n",
    "$ \\boldsymbol{y}  = X\\beta + u$ \n",
    "\n",
    "$ \\boldsymbol{u} = \\boldsymbol{y} - X\\beta $\n",
    "\n",
    "$\\boldsymbol{u}^2 = \\boldsymbol{u}\\boldsymbol{u} = (\\boldsymbol{y} - X\\beta)(\\boldsymbol{y} - X\\beta)$\n",
    "\n",
    "$\\boldsymbol{u}\\boldsymbol{u} = \\boldsymbol{y}\\boldsymbol{y} - 2\\boldsymbol{y}X\\beta + X\\beta X\\beta$\n",
    "\n",
    "Computing the partial derivative of the squared error with respect to $\\beta$, and setting it to zero yields:\n",
    "\n",
    "$\\frac{ \\partial(\\boldsymbol{u}\\boldsymbol{u})}{ \\partial\\beta} = 0$\n",
    "$ = -2X\\boldsymbol{y} + 2X^2\\beta$\n",
    "\n",
    "Solving for $\\beta$ yields:\n",
    "\n",
    "$\\beta = \\frac{XY}{X^2}$\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "    We do <b>not</b> reduce the $X$ values - they are matrices, and the equations here will change slightly in the next notebook when we examine these values as matrices\n",
    "</div>\n",
    "\n",
    "Below, I take you through the code to do the same. There is a bit of matrix manipulation happening here that may not be familiar, however we will cover that in the next notebook. Ignoring the transposing (`.T`), inverting (`np.linalg.inv`) and some rearranging if you are not familiar with matrix maths, this is the same result as before. Also note that generally you do *not* do this yourself, instead letting an established library do it for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_prices = prices['MedianPrice'].values\n",
    "cash_rate = prices['CashRate'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(53,)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cash_rate.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(53,)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ones = np.ones((len(cash_rate)))\n",
    "ones.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(53, 2)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.vstack([ones, cash_rate]).T\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(53, 2)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = house_prices.reshape((len(house_prices), 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((53, 2), (53, 1))"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape, Y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[975.51254806],\n",
       "       [-78.60228844]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(Y))\n",
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those are the new `m` and `c` values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1d3a72b4a50>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdCUlEQVR4nO3de3zO9f/H8ce1sVnYhrI5TE6VlFNOzaEDiyR9i28HEUXpMEW+nVQ6Ccm39FUi6oty6mAkhRzmVGuYMJQoImyKtply2j6/P97fXT9jY9d2Xftch+f9dvvc5vp8Ptfnen0U18v783693g7LsixEREREAliQ3QGIiIiI2E0JkYiIiAQ8JUQiIiIS8JQQiYiISMBTQiQiIiIBTwmRiIiIBDwlRCIiIhLwytgdgC/Izc1l//79VKxYEYfDYXc4IiIiUgSWZXHkyBGqV69OUNC5x4CUEBXB/v37iYmJsTsMERERKYa9e/dSs2bNc56jhKgIKlasCJjf0PDwcJujERERkaLIysoiJibG+T1+LkqIiiDvMVl4eLgSIhERER9TlOkumlQtIiIiAU8JkYiIiAQ8JUQiIiIS8JQQiYiISMBTQiQiIiIBTwmRiIiIBDwlRCIiIhLwlBCJiIhIwFNCJCIiIgFPnap9UU4OrF4NBw5AtWrQvj0EB9sdlYiIiM9SQuRrEhJg0CD47bf/31ezJvznP9C9u31xiYiI+DA9MvMlCQnwz3/mT4YA9u0z+xMS7IlLRETExykh8hU5OWZkyLLOPpa3b/Bgc56IiIi4RAmRr1i9+uyRodNZFuzda87zlJwcWLECZs0yP5V8iYiIn9AcIl9x4IB7z3OV5i6JiIgf0wiRr6hWzb3nuUJzl0RExM8pIfIV7dubERmHo+DjDgfExJjz3Elzl0REJAAoIfIVwcHm8RScnRTlvX7rLff3I3LH3CXNPRIRES+nhMiXdO8On30GNWrk31+zptnvibk8JZ27lJAAtWvD9dfD3Xebn7Vr6zGbiIh4FU2q9jXdu8M//lF6napLMncpb+7RmY/b8uYeeSqJExERcZGtI0SrVq2iW7duVK9eHYfDwbx58/IdtyyLF154gWrVqhEWFkZcXBw7duzId87hw4fp1asX4eHhREZG0r9/f7Kzs/Ods3nzZtq3b0+5cuWIiYnh9ddf9/SteVZwMFx3HfTsaX56ctmO4s5d0twjERHxIbYmREePHqVJkyaMHz++wOOvv/4648aNY+LEiSQnJ1O+fHk6d+7MsWPHnOf06tWLrVu3smTJEhYsWMCqVasYMGCA83hWVhadOnXi4osvJiUlhTFjxvDSSy8xadIkj9+fXyju3CV39U3S/CMRESkNlpcArLlz5zpf5+bmWtHR0daYMWOc+zIyMqzQ0FBr1qxZlmVZ1rZt2yzAWrdunfOchQsXWg6Hw9q3b59lWZb17rvvWpUqVbKOHz/uPOfpp5+2LrvsskJjOXbsmJWZmenc9u7dawFWZmamu27X98yZY1k1a1qWSWXMFhNj9hdk5sz85xa2zZzp2mfWrFn4Z4qIiJwmMzOzyN/fXjupeteuXaSlpREXF+fcFxERQevWrUlKSgIgKSmJyMhIWrRo4TwnLi6OoKAgkpOTnedcc801hISEOM/p3Lkz27dv588//yzws0eNGkVERIRzi4mJ8cQt+pbu3WH3bkhMhJkzzc9duwqfA1TSvknqfSQiIqXIaxOitLQ0AKKiovLtj4qKch5LS0ujatWq+Y6XKVOGypUr5zunoGuc/hlnGjp0KJmZmc5t7969Jb8hf+DK3KWS9E3S/CMRESllXpsQ2Sk0NJTw8PB8m7ioJH2T1PtIRERKmdcmRNHR0QCkp6fn25+enu48Fh0dzcGDB/MdP3XqFIcPH853TkHXOP0zxEOK2zdJvY9ERKSUeW1CVKdOHaKjo1m2bJlzX1ZWFsnJycTGxgIQGxtLRkYGKSkpznOWL19Obm4urVu3dp6zatUqTp486TxnyZIlXHbZZVSqVKmU7iaAuTr3CNzT+0hzj0RExAUOyypookbpyM7OZufOnQA0a9aMN998k+uvv57KlStTq1YtRo8ezWuvvca0adOoU6cOw4YNY/PmzWzbto1y5coB0KVLF9LT05k4cSInT57kvvvuo0WLFsycOROAzMxMLrvsMjp16sTTTz/Nli1b6NevH2PHjs1Xnn8uWVlZREREkJmZqcdnpSEnx4zo7NtX8Dwih8OMMu3alf+RW977CnvcVtj7RETEL7n0/e3xmrdzSExMtICztr59+1qWZUrvhw0bZkVFRVmhoaFWx44dre3bt+e7xqFDh6yePXtaFSpUsMLDw6377rvPOnLkSL5zNm3aZLVr184KDQ21atSoYb322msuxelK2Z64yZw5luVwmO30svu8fQWV3icmFq3UPzGxtO9GRERs4Mr3t60jRL5CI0Q2SUgw1Wanj/jExJjJ2AU9cps1y8wZOp+ZM02lnIiI+DVXvr+1lpl4L1fXbTujBUOhinqeiIgEDCVE4t3yeh+JiIh4kNdWmYm47IwWDCU+T0REAoYSIpvt2QOnTtkdhZ8o6XIhIiISsJQQ2Sg7Gzp0gGuvhZ9/tjsaP9CmzfnL6YODzXkiIiKnUUJko82b4fff4dtvoUkTeP/9gtvuSBF9++35l+jIyTHniYiInEYJkY3atDFJ0bXXwtGj8MADpqjqjJVGpKhKuuSHiIgELCVENrv4Yli+HP79bwgJgS++gEaNYP58uyPzQSq7FxGRYlJC5AWCguBf/4L1600y9PvvZqTo/vvhyBG7oxMREfF/Soi8SKNGsG4dPPmkWXbrgw+gaVP45hu7I/MRKrsXEZFiUkLkZUJD4fXXzaLwtWrBL7/ANdfAc8/BiRN2R+flVHYvIiLFpITIS117rZlw3bcv5ObCyJFw9dWwbZvdkXkxld2LiEgxKSHyYhERMHUqfPYZVKkC338PzZvDuHEmSZIzqOxeRESKSQmRD+jRA1JToUsXOHbMLADfuXP+ReAFld2LiEixKSHyEdWqwZdfwrvvQlgYLF1qJmHPnm13ZF5Ec4hERKSYlBD5EIcDHn4YNm6Eli0hIwN69oS774Y//7Q7Oi+gOUQiIlJMSoh80KWXmlL8l14y3++zZpnRoqVL7Y7MZppDJCIixaSEyEeVLQsvvmi+2y+5BPbtgxtugMGD4e+/7Y7OJppDJCIixaSEyMe1amWqzx5+2Lz+z3+gRQuzL+BoDpGIiBSTEiI/UL68mWz91VcQHW16FbVuDaNGnf8Jkl9p3970JziXKlXMeSIiIqdRQuRHunQx5fndu8PJk/Dss6bB4y+/2B2ZiIiId1NC5GcuvNA0cpw2DSpWNJOvmzQx66JZlt3Redjq1XDo0LnPOXTInCciInIaJUR+yOGAPn3M0h/t20N2Ntx/P9x2m5+va6pJ1SIiUkxKiPxY7dpmkdjXXzdVaZ9/bsrzv/jC7sg8RJOqRUSkmJQQ+bngYHjySVi3Dq680owQ3XILDBhgRo78Svv2ULOmGSIriMMBMTGaVC0iImdRQhQgmjQxSdETT5i8YPJks8+vehQGB5u+A4VNlrIseOut83ezFhGRgKOEKICUKwdjxsDy5Wag5JdfzGDJ88+bqjQREZFApYQoAF13nZlwfc89kJsLI0ZAbCz88IPdkZVQTg4MGlT4cYfDtPIOqOZMIiJSFEqIAlRkJHz4IXzyCVSuDCkpcNVV8PbbJknySatXw2+/FX7csmDvXpXdi4jIWZQQBbjbbzfNHDt1gmPH4LHH4MYbzdpoPkdl9yIiUkxKiITq1WHRInjnHTPPaMkSU57/8cd2R+Yild2LiEgxKSESwEyviY83i8K2aAF//gl33QW9e0NGht3RFZHK7kVEpJiUEEk+DRqYUvxhwyAoCGbMMKNFy5fbHVkR5JXdw9lJUd5rld2LiEgBlBDJWcqWhVdeMeug1atn5il37AhDhph5Rl6te3fTbCnojP+1g4LM/u7d7YlLRES8mhIiKdTVV8PGjaarNcDYseZx2saNdkZ1HgkJ8O9/n11an5Nj9ick2BOXiIh4NSVEck4VKsB775n1z6pWha1boVUrGD3aC9v55PUhKqxTNagPkYiIFEgJkRTJzTfDli1w662mq/Uzz5gGj7t22R3ZadSHSEREikkJkRTZRReZJ07//a8ZOVqzBho3hilTzj0oU2rUh0hERIpJCZG4xOGA++4zS3+0awfZ2dCvH/ToAb//bnNw6kMkIiLFpIRIiqVOHVixAl57zVSlzZ1ryvO//NLGoNSHSEREisnrE6IjR44wePBgLr74YsLCwmjTpg3r1q1zHrcsixdeeIFq1aoRFhZGXFwcO3bsyHeNw4cP06tXL8LDw4mMjKR///5kZ2eX9q34neBgePppSE6Ghg0hPd3MNXroITNyZEtA6kMkIiLF4PUJ0f3338+SJUv46KOPSE1NpVOnTsTFxbHvf4ttvf7664wbN46JEyeSnJxM+fLl6dy5M8dOa5jTq1cvtm7dypIlS1iwYAGrVq1iQF4tuZRYs2ZmcdjHHzev33vP7PvuOxuC6d4dPvsMatTIv79GDbNffYhERKQglhf766+/rODgYGvBggX59l911VXWc889Z+Xm5lrR0dHWmDFjnMcyMjKs0NBQa9asWZZlWda2bdsswFq3bp3znIULF1oOh8Pat29fkeLIzMy0ACszM9MNd+Xfli61rJo1LQssKyjIsoYNs6wTJ0o5iDlzLKtGDRNE3lajhtkvIiIBw5Xvb68eITp16hQ5OTmUK1cu3/6wsDDWrFnDrl27SEtLIy4uznksIiKC1q1bk5SUBEBSUhKRkZG0aNHCeU5cXBxBQUEkJycX+LnHjx8nKysr3yZF07EjpKZCr16QmwvDh0NsLPz4YykFkJAA//wn/G8E0Wn/frNfjRlFRKQAXp0QVaxYkdjYWIYPH87+/fvJyclh+vTpJCUlceDAAdLS0gCIiorK976oqCjnsbS0NKpWrZrveJkyZahcubLznDONGjWKiIgI5xYTE+OBu/NfkZEwfTrMnm1+nZJiHqG9846Hy/PP1Zgxb58aM4qISAG8OiEC+Oijj7Asixo1ahAaGsq4cePo2bMnQWeuVeVGQ4cOJTMz07nt3bvXY5/lz+680zRzvOEGswbao4/CjTeawRqPUGNGEREpJq9PiOrVq8fKlSvJzs5m7969rF27lpMnT1K3bl2io6MBSE9Pz/ee9PR057Ho6GgOHjyY7/ipU6c4fPiw85wzhYaGEh4enm+T4qlRAxYtgnHjoFw5+PprU57/6ace+DA1ZhQRkWLy+oQoT/ny5alWrRp//vknixcv5h//+Ad16tQhOjqaZcuWOc/LysoiOTmZ2NhYAGJjY8nIyCAlJcV5zvLly8nNzaV169alfh+BKCjIjA5t2ABXXQWHD8Mdd8A990BGhhs/SI0ZRUSkmByW5RWLLhRq8eLFWJbFZZddxs6dO3nyyScpV64cq1evpmzZsowePZrXXnuNadOmUadOHYYNG8bmzZvZtm2bczJ2ly5dSE9PZ+LEiZw8eZL77ruPFi1aMHPmzCLFkJWVRUREBJmZmRotKqETJ8xE65EjzaTrmBj48EOzLlqJ5eRA7dpmQnVB/1s7HKZx465d6kUkIhIAXPn+9voRoszMTOLj42nQoAF9+vShXbt2LF68mLJlywLw1FNP8eijjzJgwABatmxJdnY2ixYtyleZNmPGDBo0aEDHjh256aabaNeuHZMmTbLrlgJaSIhJiNasgXr1zJSeDh3giSfMPKMSUWNGEREpJq8fIfIGGiHyjOxsGDIEJk82rxs1MtVpjRuX8MIJCaba7PQJ1jVrmmRJjRlFRAKGX40Qif+qUAEmTYL586FqVdO/qGVLGDPGDZXxZ+b5yvtFROQclBCJ7bp1M8nQLbeYOUZPPWUeo+3eXYyLqTGjiIgUgxIi8QpVq8K8efDBB2bkaNUq8+jsww9dGNxRY0YRESkmJUTiNRwO6NcPNm2CNm3gyBHo2xduvx3++KMIF1BjRhERKSYlROJ16tY1I0QjR0KZMjBnjplwvXDhed6oxowiIlJMSojEKwUHw9ChkJwMl18OaWlw003wyCNw9Gghb1JjRhERKSYlROLVrrrKLA47aJB5PWGCWSg2ObmAk9u3N+X1Z/YgyuNwmE6Q7dt7LF4REfFNSojE64WFmX6KS5aYtdF27IC2beGll+DkydNOVGNGEREpJiVE4jPi4kx5fs+eplDs5ZdNYrR9+2knde8On31mMqfT1axp9helMWNODqxYAbNmmZ+qShMR8XvqVF0E6lTtfWbNMvOJMjLMCNIbb8BDD502MJSTY6rJDhwwc4baty/ayJC6XIuI+A1Xvr+VEBWBEiLv9NtvcO+9sGyZed2li+ljVOw503lNHc/8I5GXZRV1hElERLyClu6QgFCzJnz9tZkWFBpqyvIbNTJl+i5TU0cRkYCmhEh8WlCQyWM2bDDVZ4cOmUGevn0hM9OFC6mpo4hIQFNCJH6hYUP47jt49lmTJH34oVn6Y+XKIl5ATR1FRAKaEiLxGyEhMGKE6XJdty7s2QPXX28Wiz1+/DxvVlNHEZGApoRI/E7btrBxI/Tvb550jRkDLVvC5s3neJOaOoqIBDQlROKXKlaE99+HefPgootM/6KWLeHf/y5kXrSaOoqIBDQlROLX/vEPkwzdfDOcOAFPPgkdO8KvvxZwckmbOqqho4iIz1IfoiJQHyLfZ1mmR9HgwWZx2PBweOcd6N27gKdkxWnqqIaOIiJeR40Z3UwJkf/YuRP69IGkJPP6n/+EiROhSpUSXFQNHUVEvJIaM4oUon59U4X26qtQpozJVRo1gkWLinlBNXQUEfELSogk4JQpA889Z/oWNWhgnox16QLx8fDXXy5eTA0dRUT8ghIiCVjNm0NKCjz6qHn97rum2/W6dS5cRA0dRUT8ghIiCWgXXADjxsHixVC9Ovz0E8TGwiuvwKlTRbiAGjqKiPgFJUQiQKdOpjz/jjvMdJ8XX4R27WDHjvO8UQ0dRUT8ghIikf+pXBlmz4YZMyAiApKToWlTeO+9gudMA2roKCLiJ5QQiZzG4YC77zajRR06mEnWDz1kGjumpRXypsIaOtaooZJ7EREfoYRIpAAxMbBkCYwdC6Gh8NVXcOWVMHfuOd505jCSWnyJiPgMJUQihQgKMi2E1q+HJk3g0CEz2NOvH2RlnXZiXmPGffvyX2D/frM/IaE0wxYRkWJQQiRyHldeaeYTPfOMeaQ2ZYpJkFavRo0ZRUT8hBIikSIIDYVRo2DlSqhdG3bvhmuvhWfu2cfx3w4W/kY1ZhQR8QlKiERc0L49bNpkHptZFoyeVYvWJLOFK879RjVmFBHxakqIRFwUHg4ffGCmBl0YcYJNNKUF6xnLYHIppB+RGjOKiHg1JUQixXTbbZC6NZiu5ZZynHIMYSxxLGUvNf//JDVmFBHxCUqIREogukYwX0zP4j0e5AKOkkgHGpHKDO7GQo0ZRUR8hRIikRJy9OjOgDmd2RR9I1eTRCaR9GYGPcPmcXjK52rMKCLiA5QQibhD9+7U/20Fq5eeYPjtmykTnMvHf99Co2e7sWSJ3cGJiMj5KCEScaMywRbP37aVpPHfc9llFvv3m4VjH3vMLAMiIiLeSQmRiDskJJgGRddfD3ffTYuHWrDhyKUM7PIzAG+/Dc2bQ0qKvWGKiEjBlBCJlFTe0h2//ZZv9wUHfubtRZew6Pk1VKsGP/4IV18Nr74Kp07ZFKuIiBRICZFISRRh6Y7O0+4mdWMOt99uEqFhw+Caa2DnzlKOVURECuXVCVFOTg7Dhg2jTp06hIWFUa9ePYYPH4512pePZVm88MILVKtWjbCwMOLi4tixY0e+6xw+fJhevXoRHh5OZGQk/fv3Jzs7u7RvR/zR6tVnjQzl87+lO6psW83HH8NHH5nGjklJ0LQpTJ5ccC4lIiKly6sTotGjRzNhwgTeeecdfvjhB0aPHs3rr7/O22+/7Tzn9ddfZ9y4cUycOJHk5GTKly9P586dOXbsmPOcXr16sXXrVpYsWcKCBQtYtWoVAwYMsOOWxN8UdUmOAwdwOKB3b0hNheuug6NHYcAAuOUWSE/3aJQiInI+lhfr2rWr1a9fv3z7unfvbvXq1cuyLMvKzc21oqOjrTFjxjiPZ2RkWKGhodasWbMsy7Ksbdu2WYC1bt065zkLFy60HA6HtW/fviLFkZmZaQFWZmZmSW9J/E1iomWZQZ5zb4mJ+d6Wk2NZb7xhWSEh5vCFF1rW3Ll23ICIiP9y5fvbq0eI2rRpw7Jly/jpp58A2LRpE2vWrKFLly4A7Nq1i7S0NOLi4pzviYiIoHXr1iQlJQGQlJREZGQkLVq0cJ4TFxdHUFAQycnJBX7u8ePHycrKyreJFKh9e6hZ0yzRUZBClu4ICoIhQ2D9emjcGP74wywF0r8/HDlSCnGLiEg+Xp0QPfPMM9x11100aNCAsmXL0qxZMwYPHkyvXr0ASEtLAyAqKirf+6KiopzH0tLSqFq1ar7jZcqUoXLlys5zzjRq1CgiIiKcW0xMjLtvTfxFcDD85z/m12cmRY7zL93RqBGsXQtPPWVO/+9/oUkT+OYbz4UsIiJn8+qE6JNPPmHGjBnMnDmTDRs2MG3aNP79738zbdo0j37u0KFDyczMdG579+716OeJj+veHT77DGrUyL+/Zk2z/zxLd4SGwujRsGIFXHwx7NplqtCefRZOnPBc2CIi8v/K2B3AuTz55JPOUSKARo0a8euvvzJq1Cj69u1LdHQ0AOnp6VSrVs35vvT0dJo2bQpAdHQ0Bw8ezHfdU6dOcfjwYef7zxQaGkpoaKgH7kj8VvfucPPN8O678PPPUK8ePPIIhIQU+RLXXAObN5uu1tOmwahRsHAhTJ8OV1zhwdhFRMS7R4j++usvgoLyhxgcHExubi4AderUITo6mmXLljmPZ2VlkZycTGxsLACxsbFkZGSQclqL4OXLl5Obm0vr1q1L4S4kICQkmCTo8cfhnXfMz3r1zH4XhIfD1KlmYKlKFdi40XS4fust+N//9iIi4gFenRB169aNESNG8OWXX7J7927mzp3Lm2++yW233QaAw+Fg8ODBvPrqq8yfP5/U1FT69OlD9erVufXWWwG4/PLLufHGG3nggQdYu3Yt33zzDQMHDuSuu+6ievXqNt6d+I1COlWzb5/Z72JSBNCjhynP79IFjh83+VWnTudueSQiIiVQClVvxZaVlWUNGjTIqlWrllWuXDmrbt261nPPPWcdP37ceU5ubq41bNgwKyoqygoNDbU6duxobd++Pd91Dh06ZPXs2dOqUKGCFR4ebt13333WkSNHihyHyu6lUKdOWVbNmoWX2zsclhUTY84rhtxcy5owwbIuuMBcLjLSsmbOdPM9iIj4KVe+vx2WpT6555OVlUVERASZmZmEh4fbHY54kxUrzIKu55OYaLoxFtNPP8E995iKNIC77jLTlSpVKvYlRUT8nivf3179yEzE67nQqbokLr3UlOK//LKp4J8925TsL11aosuKiMj/KCESKYnTqhvdct45lCkDL7xg1kG79FIzRemGG2DwYPj77xJfXkQkoCkhEimJYnaqLomWLeH77yE+3rz+z3+gRQvYsMFtHyEiEnCUEImURAk7VRfXBReY6v6FCyE6GrZtg9atYeRIyMlx60eJiAQEJUQiJVXCTtUlceONsGWLKdM/dQqee840ePzlF499pIiIX1KVWRGoykyKJCcHVq82E6irVTOPydw8MlQYyzIdrQcOhKwsqFDBDEz161f40zwREX/n0SqzadOm8eWXXzpfP/XUU0RGRtKmTRt+/fVX16MV8RfBwaa0vmdP89OVZCgnx5Twz5plfrr43MvhMGX5mzebEaLsbLj/frjtNjhj5RoRESmAywnRyJEjCQsLAyApKYnx48fz+uuvc+GFF/L444+7PUARv5eQALVrm35Gd99tftauXawO1xdfDMuXw5gxZhm1zz835flffOH2qEVE/IrLCdHevXupX78+APPmzaNHjx4MGDCAUaNGsXr1arcHKOLXPLDsR3AwPPEErFtnkqGDB+GWW2DAADNyJCIiZ3M5IapQoQKHDh0C4Ouvv+aGG24AoFy5cvytZigiRZeTA4MGmQlAZ8rbN3hwscvGGjc2na2feMI8Ups8GZo0gW+/LX7IIiL+yuWE6IYbbuD+++/n/vvv56effuKmm24CYOvWrdSuXdvd8Yn4r9Wrz71aq2XB3r3mvGIqV848Plu+HGrVMtVn7dvD88/DiRPFvqyIiN9xOSEaP348sbGx/P7778yZM4cqVaoAkJKSQs+ePd0eoIjfKqVlP8DM8d682Uy8zs2FESMgNhZ++KHElxYR8Qsquy8Cld2LR5TSwrBn+vRTeOghOHzYjCCNHm3K9YPUlUxE/IzHF3ddvXo1vXv3pk2bNuzbtw+Ajz76iDVr1hTnciKByYZlPwBuvx1SU6FzZzh2zExj6tz53E/vRET8ncsJ0Zw5c+jcuTNhYWFs2LCB48ePA5CZmcnIkSPdHqCI37Jp2Q+A6tXNsh/vvANhYbB0qalI+/hjt3+UiIhPcDkhevXVV5k4cSKTJ0+mbNmyzv1t27Zlg1aXFHGNjct+OBxmgdgNG8zisBkZcNdd0KsX/Pmnxz5WRMQruZwQbd++nWuuueas/REREWRkZLgjJpHA0r077N5t5grNnGl+7tpVtGSohB2uARo0MKX4L7xgBqNmzjQl+8uWuXwpERGf5XJCFB0dzc6dO8/av2bNGurWreuWoEQCTnGW/XBjh+uyZeHll2HNGqhf38wniouDIUPMPCMREX/nckL0wAMPMGjQIJKTk3E4HOzfv58ZM2bwxBNP8PDDD3siRhE5kwc6XANcfTVs3Giq0ADGjoXmzeH770sWroiIt3O57N6yLEaOHMmoUaP466+/AAgNDeWJJ55g+PDhHgnSbiq7F6+Sk2NGggorC3M4zBykXbtKNCH7yy+hf39ITzcjSK+8Ak8+6ZE53iIiHuHK93ex+xCdOHGCnTt3kp2dTcOGDalQoUKxgvUFSojEq5Ri/6LffzdroM2bZ163bQsffQR16pTosiIipcKjfYgyMzM5fPgwISEhNGzYkFatWlGhQgUOHz5MVlZWsYMWkSIqxQ7XF11knr5NmQIVK8I335gJ11OmFLwEm4iIr3I5IbrrrruYPXv2Wfs/+eQT7rrrLrcEJSLnUK2ae887D4cD7r0XNm2Cdu0gOxv69TNFcL//7paPEBGxncsJUXJyMtcXMFx/3XXXkZyc7JagROQcbOpwXaeOeVr32mtmTtG8eXDllbBggVs/RkTEFi4nRMePH+fUqVNn7T958iR///23W4ISkXOwscN1cDA8/TSsXQtXXAEHD0K3bvDgg2bkSETEV7mcELVq1YpJkyadtX/ixIk0b97cLUGJyHnY2OEaoGlTWL/e9CkCmDTJ7EtK8ujHioh4jMtVZt988w1xcXG0bNmSjh07ArBs2TLWrVvH119/TXs3D9N7A1WZidc6cQLefRd+/hnq1YNHHoGQkFINYfly6NvXdAEICoJnnzVdr09b2UdExBYeL7vfuHEjY8aMYePGjYSFhdG4cWOGDh3KJZdcUuygvZkSIvFKCQlmqfrT+xHVrGkep3l4hOhMGRkwcCDMmGFeN28O06ebZUFEROxSKn2IAokSIvE6eZ2qz/zjmzeHqBQemxXk44/h4YfN4rDlysGYMWYB2cLmf4uIeJLbE6KsrCznhc7Xa8gfEwYlROJVSqlTdXHt2wf33QdLlpjXnTqZvkXVq5d6KCIS4NzemLFSpUocPHgQgMjISCpVqnTWlrdfRDxs9erCkyEwo0Z795rzbFCjBixaBG+/bUaJvv7alOd/+qkt4YiIFEmZopy0fPlyKleuDEBiYqJHAxKR8yjFTtXFFRRk5hTFxUHv3pCSAnfcYX799tsQGWlbaCIiBSpSQnTttdcCcOrUKVauXEm/fv2oWbOmRwMTkUKUcqfqkmjQwJTiDx8OI0aYidYrV8K0aUVbjk1EpLS41IeoTJkyjBkzpsDGjCJSSmzqVF1cZcvCK6/AmjWmM8DevdChA/zrX3DsmN3RiYgYLjdm7NChAytXrvRELCJSFDZ2qi6J2FjYuBEGDDCv33wTWrY0a6SJiNitSI/MTtelSxeeeeYZUlNTad68OeXLl893/JZbbnFbcCJSiO7d4YknTFaRk/P/+4OCTPtoG0rui6JCBXjvPbPcR//+sGWLSYpefdWMGHlZDiciAcTlPkRBQYUPKjkcDnJO/8vZT6jsXrxOYX2IwIwS2dSHyBW//w4PPACff25et28PH35oOgqIiLiD28vuT5ebm1vo5o/JkIjXyckxHarP9W+ZwYPzjxx5oYsugrlz4YMPzMjR6tXQuDFMnXruWxMR8QSXEqLdu3czefJk3n33XbZu3eqpmETkXLy8D5ErHA7o18/MI2rbFo4cMU0de/SAP/6wOzoRCSRFTogSExO54oorePDBBxk4cCDNmjVj+vTpnoxNRAriA32IXFW3rinHHzXKVKXNnWuaOX71ld2RiUigKHJCNGzYMG644Qb27dvHoUOHeOCBB3jqqac8GRsAtWvXxuFwnLXFx8cDcOzYMeLj46lSpQoVKlSgR48epKen57vGnj176Nq1KxdccAFVq1blySefVOsA8V0+1IfIFcHB8MwzkJwMDRtCejp07WrWRjt61O7oRMTfFXlSdWRkJN9++y0NGzYE4K+//iI8PJz09HSqVKnisQB///33fHOTtmzZwg033EBiYiLXXXcdDz/8MF9++SVTp04lIiKCgQMHEhQUxDfffANATk4OTZs2JTo6mjFjxnDgwAH69OnDAw88wMiRI4sUgyZVi1fJW8ts377CJ1XbuJaZO/z9Nzz7rOkeAHDJJfDRR9C6ta1hiYiPcen72yoih8Nhpaen59tXoUIF6+effy7qJdxi0KBBVr169azc3FwrIyPDKlu2rPXpp586j//www8WYCUlJVmWZVlfffWVFRQUZKWlpTnPmTBhghUeHm4dP368SJ+ZmZlpAVZmZqZ7b0akuObMsSyHw2wmLTJb3r45c+yO0C2WLrWsGjXMrQUHW9YLL1jWiRN2RyUivsKV72+XJlUvXryY+fPnO7fc3FyWLVuWb58nnThxgunTp9OvXz8cDgcpKSmcPHmSuLg45zkNGjSgVq1aJCUlAZCUlESjRo2IiopyntO5c2eysrIKnRh+/PhxsrKy8m0iXqV7d1NaX6NG/v01avhEyX1RdewIqanQs6cZGHvlFTP5evt2uyMTEX/jUmPGvn37nrXvwQcfdP7a032I5s2bR0ZGBvfeey8AaWlphISEEHnGSpFRUVGkpaU5zzk9Gco7nnesIKNGjeLll192b/AinnDmIzM/rFevVAlmzjTNHB95BNatg2bN4N//NvOLClvBRETEFUUeITpX/6HS6kP0wQcf0KVLF6pXr+7Rzxk6dCiZmZnObe/evR79PBGX5TVm3Lcv//79+83+hAR74vKgnj3NaFHHjmaOUXw83HSTTxXTiYgXc7kxo11+/fVXli5dyv333+/cFx0dzYkTJ8jIyMh3bnp6OtHR0c5zzqw6y3udd86ZQkNDCQ8Pz7eJeI1zNWbM2+cDjRmLo2ZN+Pprs5RbuXKwaJEpz58zx+7IRMTX+UxCNGXKFKpWrUrXrl2d+5o3b07ZsmVZtmyZc9/27dvZs2cPsbGxAMTGxpKamsrBgwed5yxZsoTw8HBnxZyIT/GjxozFERQEjz0GKSnm0dnhw2ZQrG9fyMy0OzoR8VU+kRDl5uYyZcoU+vbtS5ky/z/tKSIigv79+zNkyBASExNJSUnhvvvuIzY2lquvvhqATp060bBhQ+655x42bdrE4sWLef7554mPjyc0NNSuWxIpPj9szFgcDRvCd9/Bc8+ZJOnDD83SHytX2h2ZiPgin0iIli5dyp49e+jXr99Zx8aOHcvNN99Mjx49uOaaa4iOjibhtPkTwcHBLFiwgODgYGJjY+nduzd9+vThlVdeKc1bEHEfP23MWBwhIfDqq7Bqlel2vWcPXH89PPkkHD9ud3Qi4ktcXu0+EKkxo3iVAGjMWBxHjsCQIfD+++Z1o0YwfboZNRKRwOTR1e4BMjIyeP/99xk6dCiHDx8GYMOGDew7s+JFRNwvONjMKi7s3zKWZVo8B1AyBFCxIkyeDJ9/DhddZCrSWraEMWP8cn65iLiZywnR5s2bufTSSxk9ejT//ve/nRVeCQkJDB061N3xiYi45JZbYMsW07foxAl46ino0AF277Y7MhHxZi4nREOGDOHee+9lx44dlCtXzrn/pptuYtWqVW4NTkQKkFd2XxiHw2/L7ouqalUzUjR5MpQvb+YYNW4M06b5Ze9KEXEDlxOidevW5etOnadGjRqFdn4WETcK8LL7onI44P77YdMmaNPGzDG69164/Xb44w+7oxMRb+NyQhQaGlrg2l4//fQTF110kVuCEpFzUNm9S+rVM6X4I0ZAmTKmiWOjRqapo4hIHpcToltuuYVXXnmFkydPAmb9sj179vD000/To0cPtwcoImdQ2b3LypSBZ5+F5GS4/HJIS4MuXczyH0eP2h2diHgDlxOiN954g+zsbKpWrcrff//NtddeS/369alYsSIjRozwRIwicrr27U1ZfWGrmjocEBNjzpN8rrrKdLjOm4L17rtm39q19sYlIvYrdh+ib775hk2bNpGdnc1VV11FXFycu2PzGupDJF4nIQHONSI7Zw5071568figpUvNnKJ9+0yHgmHDTNfr05rhi4iPc+X7u9h/9Nu2bUvbtm0BzlpcVUTE28XFwebN8Mgj8PHH8NJL8NVX8NFHcOmldkcnIqXN5Udmo0eP5uOPP3a+vuOOO6hSpQo1atRg06ZNbg1ORAqgsnu3qVwZZs+GmTMhIsI8OmvWDCZOVHm+SKBxOSGaOHEiMTExgFk1fsmSJSxcuJAuXbrw5JNPuj1AETmDyu7drmdP09m6Qwf46y94+GHo2tVMvhaRwOByQpSWluZMiBYsWMAdd9xBp06deOqpp1i3bp3bAxSRM6js3iNiYmDJEhg7FkJDYeFCuPJKM11LRPyfywlRpUqV2Lt3LwCLFi1yTqa2LIscDdGLeJ67yu5zcmDFCpg1y/zUn1+CgszTxpQUaNoUDh0yc9fvuw8KaL8mIn7E5YSoe/fu3H333dxwww0cOnSILl26APD9999Tv359twcoImdwR9l9QgLUrg3XXw93321+1q6t4ZD/ueIK07No6FDz2zl1qln6Q6sTifgvlxOisWPHMnDgQBo2bMiSJUuoUKECAAcOHOCRRx5xe4Aicoa81e7h7KQo7/W5VrtPSIB//vPseUj79pn9SooACAmBkSNNElS7Nvz6K1x3HTz9NBw/bnd0IuJuxe5DFEjUh0i8UkKCqTY7PbGJiTHJUGE9iHJyzLd7YZOyHQ4z+rRrV+EJVQDKyoLHH4f//te8btIEpk83c4xExHu58v3tckL04YcfnvN4nz59XLmcT1BCJF4rJ8dUkx04YOYMtW9/7kRmxQrzeOx8EhPNcIjkM28ePPCAWRw2JARGjTJzjoJcHmsXkdLg0YSoUqVK+V6fPHmSv/76i5CQEC644AIOHz7sesReTgmR+I1Zs8ycofOZOdPUohfE1STMz6Snw/33w4IF5vX115s5RrVq2RqWiBTAle9vl/9d8+eff+bbsrOz2b59O+3atWPWrFnFDlpESkFJK9Q0GZuoKJg/HyZNgvLlzWBao0bmEZomIIj4LrfNIVq/fj29e/fmxx9/dMflvIpGiMRv5M0h2rev4G/vc80hypuMfeb78iZyf/ZZwK2ftnMn3HMPfPedeX377abLdeXK9sYlIoZHR4gKU6ZMGfbv3++uy4mIJxS3Qi1vuZCCkqi8fQG4XEj9+ubp4fDhZlHYTz81o0Vff213ZCLiKpdHiObPn5/vtWVZHDhwgHfeeYeYmBgWLlzo1gC9gUaIxO+4WqGmydjntX499O4N27eb1wMHwujRcMEF9sYlEsg8utr9rbfemu+1w+HgoosuokOHDrzxxhuuXk5E7NC9O/zjH0WfHK3lQs6rRQvYsAGeeQbefhveeccsBTJ9ujkmIt7N5YQoNzfXE3GISGkLDi76aI67lgvxcxdcAOPGwc03m+U+tm+H2Fh48UWTKJVx+W9cESktJZpDZFkW6usoEgDat4cqVc59TpUq514uJIB06gSpqXDHHXDqFAwbZn5rdu60OzIRKUyxEqIPP/yQRo0aERYWRlhYGI0bN+ajjz5yd2wiIj6rcmWYPds8MouIMJVoTZqYcn39O1LE+7icEL355ps8/PDD3HTTTXzyySd88skn3HjjjTz00EOMHTvWEzGKiN1WrzZLv5/LoUPmPE/JyTGTu2fNMj99oKLN4YBevWDzZjMn/a+/4MEHoVs3SEuzOzoROZ3LT7TffvttJkyYkG+JjltuuYUrrriCl156iccff9ytAYqIF7B7UnVBVXE1a5oWAj7Q+6hWLVi61BTxDR0KX35pyvMnT4Yz6lRExCYujxAdOHCANm3anLW/TZs2HAjgChMRv2bnpOq8hpBnLki7b5/Z7yNdsoOCYMgQSEkxj87++ANuuw369TOLx4qIvVxOiOrXr88nn3xy1v6PP/6YSy65xC1BiYiXad/ejMic2cwxj8Nh+hi5e1K1HzaEvPJKSE6Gp582v21TppgEyZNPG0Xk/Fx+ZPbyyy9z5513smrVKtq2bQvAN998w7JlywpMlETED+R1uP7nP823+OkJyrk6XJfU6tVnjwydzrJg715zng81hAwNhddeg65dzdIfu3fDtdeaJOnllyEkxO4IRQKPyyNEPXr0IDk5mQsvvJB58+Yxb948LrzwQtauXcttt93miRhFxBt0727WK6tRI//+mjU9t46Z3XOXPKx9ezPh+t57TW732mvQqhVs3Wp3ZCKBx22Lu/ozLd0hcpqcnKJ3uC6pAFoyJCEBBgwwxXqhoTBqlHlaGOS2FSdFAo8r399FToiyijjrzx8TBiVEIjY5ccK0fz7XHKHgYFPP7k3PmYqZNKalQf/+8NVX5nWHDjB1qpmeJSKu88hq95GRkVSqVKnQLe+4iIjbfPvt+SdM5+SY87xFQgLUrm1Gtu6+2/ysXbtI1XDR0bBgAUyYYPLA5ctNef7MmWrmKOJpRZ5UnZiY6Py1ZVncdNNNvP/++9Q4cz6BiIi7+NocorwWAWdmL3ktAoow18rhgIcego4dzYTr5GTT3HH+fHj3XdMBW0Tcr9hziCpWrMimTZuoW7euu2PyOnpkJmITX5pDlJNjRoIKq4pzOMwE9F27ijzn6tQpM5fo5ZfN5WvUMGX6N9zgvrBF/JlHHpmJiJQ6u/ofFYcrLQKKqEwZszBsUhJceqkZaOrUyUy2/vtvN8QsIk5KiETEe+X1P4KzkyJP9j8qDg8+3mvZEr7/HuLjzetx4+Cqq0zXaxFxjxIlRI7C/tUmIuIudvQ/Kg4PL29ywQXwzjuwcKG5xI8/wtVXw4gR5tGaiJRMkecQdT/jL50vvviCDh06UL58+Xz7E3xkXSFXaA6RiBcozf5HxVGKLQIOHYIHH4Q5c8zrNm3gwyk51Nvvxb8/IjbwyByiiIiIfFvv3r2pXr36Wfvdbd++ffTu3ZsqVaoQFhZGo0aNWL9+vfO4ZVm88MILVKtWjbCwMOLi4tixY0e+axw+fJhevXoRHh5OZGQk/fv3Jzs72+2xiogHBQebidM9e5qf3vZlX4otAqpUgU8/hQ8/hPBwc8kmDY7x/vXTsVws9RcRo8hl91OmTPFkHAX6888/adu2Lddffz0LFy7koosuYseOHfn6Hb3++uuMGzeOadOmUadOHYYNG0bnzp3Ztm0b5cqVA6BXr14cOHCAJUuWcPLkSe677z4GDBjAzJkzS/2eRMRPlXKLAIfDlOVfc3QhfR8OY6V1HQ/wPvO5hfe5n6oulPqLCGB5saefftpq165docdzc3Ot6Ohoa8yYMc59GRkZVmhoqDVr1izLsixr27ZtFmCtW7fOec7ChQsth8Nh7du3r0hxZGZmWoCVmZlZzDsREb+XmGhZppbs3Ftiovs+89Qpy6pZ08rBYf2bIVYIxyywrItItz6nm2U5HJYVE2POEwlArnx/e3WV2fz582nRogW33347VatWpVmzZkyePNl5fNeuXaSlpREXF+fcFxERQevWrUlKSgIgKSmJyMhIWrRo4TwnLi6OoKAgkpOTC/zc48ePk5WVlW8TETknO1oE/K/UPwiLf/Em62lBYzbxO1X5B/N5wHqPI3v/dKnUXyRQeXVC9MsvvzBhwgQuueQSFi9ezMMPP8xjjz3GtGnTAEhLSwMgKioq3/uioqKcx9LS0qhatWq+42XKlKFy5crOc840atSofPOiYrSQkIicjx0tAs54/NaILaylFU8xGge5vM8DNGUj36486b7PFPFTXp0Q5ebmctVVVzFy5EiaNWvGgAEDeOCBB5g4caJHP3fo0KFkZmY6t71793r080TET5R2i4Az/rEHEMoJRvMMiVxPLX7lF+rR/pU4nnvOFMKJSMG8OiGqVq0aDRs2zLfv8ssvZ8+ePQBER0cDkJ6enu+c9PR057Ho6GgOHjyY7/ipU6c4fPiw85wzhYaGEh4enm8TESmS7t1h926znMjMmebnrl2lPrH5Wlaxmcb0ZSq5uQ5GjjR9i7ZtK9UwRHyGVydEbdu2Zfv27fn2/fTTT1x88cUA1KlTh+joaJYtW+Y8npWVRXJyMrGxsQDExsaSkZFBymktXZcvX05ubi6tW7cuhbsQkYBTWi0CzvjH3pkiyGIq9/HZ4NVUqWK6XTdvbjpd5+Z6JiQRX+XVCdHjjz/Od999x8iRI9m5cyczZ85k0qRJxP+vf73D4WDw4MG8+uqrzJ8/n9TUVPr06UP16tW59dZbATOidOONN/LAAw+wdu1avvnmGwYOHMhdd91F9erVbbw7EZESKuCRWUF63HyC1FTo0gWOHTNroXXufO6l10QCjVcnRC1btmTu3LnMmjWLK6+8kuHDh/PWW2/Rq1cv5zlPPfUUjz76KAMGDKBly5ZkZ2ezaNEiZw8igBkzZtCgQQM6duzITTfdRLt27Zg0aZIdtyQiYotq1eDLL+HddyEsDJYuhUaNYPZsuyMT8Q5FXrojkGnpDhHxSrNmwd13n/+8mTPN47v/2b7dNHVct8687tkTxo+H03reivgFjyzdISIipSAnB1asMMnOihXnXg6kiI/Mzjzvssvgm2/gxRfN9KZZs8xo0dKlxY5axOcpIRIR8RYJCWYNsuuvNyM/HlyTrGxZeOklkxhdcgns2wc33ACDB8Pff7v940S8nhIiERFvkJBg1h47c6Zz3ppkBSVF56kyK8p5rVub6rOHHjKv//MfaNHC7BMJJEqIRETslpNjSr8KmtKZt2/w4LMfn1WrVrTrn+e88uVhwgQz6ToqyvQqat0aRo069xM7EX+ihEhExG7/W5OsUJYFe/eevSaZm9dPu+km2LIFbrsNTp6EZ5+Fa6+FX34p4n2I73NlDpufUUIkIlKY0vpyOGNNsiKf54H10y68EObMgalToWJFM8eoSRP4738LHsASP1KKc9i8kRIiEZGClOaXQzGrxQCPrJ/mcEDfvrB5sxlcys6G/v3NyFFRpy2JjylsDttvvxU+h83PKCESETlTcSY428lD66fVrm0u9frrpirt889Nef4XX7glavEW55rDBmZ/QXPY/IwSIhHxb64+9iruBOeScEO1mKfWTwsOhiefNE0cr7zShHDLLTBggBk5Ej9wvjlsUPAcNj+jhEhE/FdxHnsVd4JzSbipWsyTmjQxSdETT5hHapMnQ9OmkJRkW0jiLvv2ufc8H6WESERc4ytVKMV97FXcCc4l4eZqMU8pVw7GjIHly6FWLfj5Z2jXDoYNM1Vp4iVc/TP6++9Fu25Rz/NRSohEpOh8pQqlJI+97Bit8UC1mCddd52ZcH3PPZCbC6++CrGx8OOPdkcmxfozetFFRbt2Uc/zUUqIRKRo3DHRuLRGl0ry2MsdozXFuU8PVIt5UkQEfPghfPIJVK4MKSnQrBm8845JksQGxf0zeub/c4Up6nm+ypLzyszMtAArMzPT7lBE7HHqlGXVrGlZJpUoeIuJMecVZs6cs69Rs6bZ724zZ5471rxt5szCY3U4zHb6+Xn7zhVzSe/z1CnLSkw0sSUmnvv31Evs22dZnTv//+126mRZv/1md1QB5nx/Rh2Owv+MuuPPt5dy5ftbI0Qicn4lrUIp7TL2kj72Ku5ojTvu00PVYp5UvTosXGhGh8LC4OuvTXn+J5/YHVkAKcmoaN4jW4ej4Ee2DodXPbL1FCVEInJ+JalCsaOM3R2PvVzt7WPHfXoRhwPi42HDBrM47J9/wp13Qu/ekJFhd3Q+ypVHryUtBvCxR7aeoIRIRM6vJFUodpSxu2uSsiujNXbcpxdq0AC+/dZUngUFwYwZZrRo+XK7I/Mxrk6OdkcxgIcafPoKJUQicn4lqUKxo4wdSv9fvHbdpxcqWxZeecWsg1a/vskTO3aEIUPg2DG7o/MBxXn06q7WDT74yNZdlBCJyPmVpArFzqaDpfkvXh9orljarr4avv8eHnzQvB471jxO27jR1rC8W3EfvfpY6wZv5LAsrV98PllZWURERJCZmUl4eLjd4YiUvpwcM1x/rkdCMTEm2TjzL9y89+7bV/Bf8g6H+ZdtQe/1JYFyn8X05ZfQr59Z+qNsWRg+3HS9DsDfinNbscI8HjufxEQzgnOmhASTUJ3+ZzUmxiRDAfLo63SufH9rhEhEzu/0KpSCnKsKJVD+5Roo91lMXbvCli1w662mq/Uzz5jv81277I7My7hjcnQAzwMqCSVEIlI0eXNyatbMvz8m5vxzcgKlgiVQ7rOYLrrIDGD8979QoQKsWQONG8OUKYUvtB5w3PHoNYDnAZWEHpkVgR6ZiZwmJ8dUSh04YP5Sbt++6H/hluS9viRQ7rMEdu2CPn1MUgRw223w3nt+vzrE+enRq1u58v2thKgIlBCJiLhfTo5ZLPaFF8xjtKgo+OAD83gtoOVVmUH+pCjv0atGG4tMc4hERMTrBQebuUTJydCwIaSnw803w0MPQXa23dGdobTW4QM9erWJRoiKQCNEIiKedewYPPusKc0H07/oo49M6b7tCqrcqlnTTKL3ZHKiR68lpkdmbqaESESkdCxbBvfea3KP4GB47jl4/nlTqm+LvMdXZ35V6vGVT9AjMxER8UkdO0JqKvTqZQZIXnkF2rSB7dttCCbA16cLNEqIRETEq0RGwvTpMHu2+fX69dCsGYwfX8rl+VqfLqAoIRIREa90551mtCguDv7+GwYOhC5dYP/+UgpA69MFFCVEIiLitWrWhMWLYdw4KFfO/LpRIzN1x+O0Pl1AUUIkIiJeLSgIHn0UUlLgqqvg8GG4/XbT2DEz04MfnLeC/LkUZQV58QlKiERExCc0bAhJSaY8PyjIlOU3bmzaAnlEcLBZ/uJc7rpLpfB+QgmRiIj4jJAQGDHCzGOuWxf27IEOHeCJJ0wvI7fKyTGNGM9l9uxzV5mVZkNHKRElRCIi4nPatIFNm+CBB0yx1xtvQKtWsHmzGz/kfFVmcO4qs4QEsy7Z9dfD3Xebn7Vrm/3idZQQiYiIT6pQASZNgvnzoWpVU5HWsqVZH80tAzElqTLLa+h4ZkK1b5/Zr6TI6yghEhERn9atm0mGbrkFTpyAp54yj9F27y7hhYtbZaaGjj5JCZGIiPi8qlVh3jx4/30oXx5WrTITrj/8sATNHPOqzPKW6TiTw1FwlZkaOvokJUQiIuIXHA7o39/MLWrTBo4cgb59TYn+H38U44LBwWYB17yLn/lhAG+9dXaVmRo6+iQlRCIi4lfq1TMjRCNHQpkyMGeOaea4cGExLta9u+kCWaNG/v01axa+sKsaOvokr06IXnrpJRwOR76tQYMGzuPHjh0jPj6eKlWqUKFCBXr06EF6enq+a+zZs4euXbtywQUXULVqVZ588klOnTpV2rciIiKlKDgYhg6F5GS4/HJIS4ObboJHHoGjR128WPfuZkJSYiLMnGl+7tpV+Cr3xX3UJrby6oQI4IorruDAgQPObc2aNc5jjz/+OF988QWffvopK1euZP/+/XQ/7X/QnJwcunbtyokTJ/j222+ZNm0aU6dO5YUXXrDjVkREpJRddZXpcD1okHk9YYJZKDY52cULBQfDddeZRo3XXXfuZozFfdQm9rK82Isvvmg1adKkwGMZGRlW2bJlrU8//dS574cffrAAKykpybIsy/rqq6+soKAgKy0tzXnOhAkTrPDwcOv48eNFjiMzM9MCrMzMzOLdiIiI2G7JEsuqUcOywLKCgy3rxRct68QJD37gnDmWVbOm+cC8LSbG7JdS4cr3t9ePEO3YsYPq1atTt25devXqxZ49ewBISUnh5MmTxMXFOc9t0KABtWrVIikpCYCkpCQaNWpEVFSU85zOnTuTlZXF1q1bC/3M48ePk5WVlW8TERHfFhdnyvN79jQV7y+/DO3awU8/eegDXX3UJrby6oSodevWTJ06lUWLFjFhwgR27dpF+/btOXLkCGlpaYSEhBAZGZnvPVFRUaSlpQGQlpaWLxnKO553rDCjRo0iIiLCucXExLj3xkRExBaVKpncZOZMiIyEtWuhaVPzKK3Y5fnn4sqjNrGVVydEXbp04fbbb6dx48Z07tyZr776ioyMDD755BOPfu7QoUPJzMx0bnv37vXo54mISOnq2dOMFnXsCH//bSZbd+2qSvhA5tUJ0ZkiIyO59NJL2blzJ9HR0Zw4cYKMjIx856SnpxMdHQ1AdHT0WVVnea/zzilIaGgo4eHh+TYREfEvNWvC11+b+c2hoaYsv1EjU6YvgcenEqLs7Gx+/vlnqlWrRvPmzSlbtizLli1zHt++fTt79uwhNjYWgNjYWFJTUzl48KDznCVLlhAeHk7Dhg1LPX4REfEuQUGmAm3DBlN9duiQWWqsb1/IzLQ7OilNXp0QPfHEE6xcuZLdu3fz7bffcttttxEcHEzPnj2JiIigf//+DBkyhMTERFJSUrjvvvuIjY3l6quvBqBTp040bNiQe+65h02bNrF48WKef/554uPjCQ0NtfnuRETEWzRsCN99Z3oXBQWZJT+aNDENHiUweHVC9Ntvv9GzZ08uu+wy7rjjDqpUqcJ3333HRRddBMDYsWO5+eab6dGjB9dccw3R0dEknLaCcHBwMAsWLCA4OJjY2Fh69+5Nnz59eOWVV+y6JRER8VIhIaa79cqVUKcO/PqrmQf91FNw/Ljd0YmnOSzLI/Pq/UpWVhYRERFkZmZqPpGISAA4cgQefxw++MC8btwYpk83c4zEd7jy/e3VI0QiIiJ2qFgR3n8f5s6FCy+EzZuhRQt44w3IzbU7OvEEJUQiIiKFuPVW2LIFbr4ZTpyAJ54wpfq//mp3ZOJuSohERETOISoK5s+HSZOgfHlYscI8QvvoIw81cxRbKCESERE5D4cDHngANm6Eq6+GrCzo0wfuuMOU6ovvU0IkIiJSRPXrw+rV8OqrUKYMfPaZmWi9aJHdkUlJKSESERFxQZky8Nxzpm9RgwZmuY8uXWDgQPjrL7ujk+JSQiQiIlIMzZtDSgo8+qh5PX686Xa9bp29cUnxKCESEREppgsugHHjYPFiqF4dfvoJYmPhlVfg1Cm7oxNXKCESEREpoU6dIDXVTLLOyYEXX4R27WDHDrsjk6JSQiQiIuIGlSvD7NkwYwZEREByMjRtCu+9p/J8X6CESERExE0cDrj7bjNa1KGDmWT90EOmsWNamt3RybkoIRIREXGzmBhYsgTGjoXQUPjqK7jySjht/XHxMkqIREREPCAoCAYPhvXroUkT08CxRw+47z7T2FG8ixIiERERD7rySjOf6JlnzCO1qVNNgrR6td2RyemUEImIiHhYaCiMGgUrV0Lt2rB7N1x7rUmSjh+3OzoBJUQiIiKlpn172LQJ+vUzlWejR0Pr1rBli92RiRIiERGRUhQeDh98YCZYX3ihSZCaN4c334TcXLujC1xKiERERGxw222mPL9rVzhxAv71L4iLgz177I4sMCkhEhERsUl0NHzxBUycaJYBSUyExo1Nc0c1cyxdSohERERs5HDAgw/Cxo1mPlFmJvTuDXfdBYcP2x1d4FBCJCIi4gUuuQTWrDELwwYHwyefQKNG8PXXdkcWGJQQiYiIeIkyZWDYMEhKgssug/37oXNneOwxswyIeI4SIhERES/TsiVs2ADx8eb122+bSrSUFHvj8mdKiERERLzQBRfAO+/AwoVQrRr8+CNcfTW8+iqcOmV3dP5HCZGIiIgXu/FGU57/z3+aRGjYMLjmGvj5Z7sj8y9KiERERLxclSpmkvVHH5nGjklJZj20yZNVnu8uSohERER8gMNhyvFTU+G66+DoURgwAG65BdLT7Y7O9ykhEhER8SG1asGyZfDGGxASAgsWmPL8zz+3OzLfpoRIRETExwQFwZAhsH696Wz9++9w661w//1w5Ijd0fkmJUQiIiI+qlEjWLsWnnrKPFL74AMzt+ibb+yOzPcoIRIREfFhoaEwejSsWAEXXwy7dpkqtGefNYvGStEoIRIREfED11wDmzdD376QmwujRpm+Rdu22R2Zb1BCJCIi4ifCw2HqVJgzx5Tqf/89XHUV/Oc/JkmSwikhEhER8TPdu5vy/C5d4PhxGDwYOnWC336zOzLvpYRIRETED1WrBl9+CRMmQFiYKdVv1AhmzbI7Mu+khEhERMRPORzw0EOwcaNZMDYjA+6+G3r2hD//tDs676KESERExM9deqkpxX/pJQgOhtmzzWjR0qV2R+Y9lBCJiIgEgLJl4cUX4dtv4ZJLYN8+uOEGGDQI/v7b7ujsp4RIREQkgLRqZarPHnnEvB43Dpo3hw0b7I3LbkqIREREAkz58jB+PHz1FURHww8/QOvWMHIknDpld3T2UEIkIiISoLp0MeX53bubROi55+Daa+Hnn+2OrPT5VEL02muv4XA4GDx4sHPfsWPHiI+Pp0qVKlSoUIEePXqQnp6e73179uyha9euXHDBBVStWpUnn3ySU4GaAouIiJzmwgvhs89g2jSoWNHMMWrSBN5/HyzL7uhKj88kROvWreO9996jcePG+fY//vjjfPHFF3z66aesXLmS/fv30717d+fxnJwcunbtyokTJ/j222+ZNm0aU6dO5YUXXijtWxAREfFKDgf06WOW/rjmGjh6FB54AG69FQ4etDu60uETCVF2dja9evVi8uTJVKpUybk/MzOTDz74gDfffJMOHTrQvHlzpkyZwrfffst3330HwNdff822bduYPn06TZs2pUuXLgwfPpzx48dzQqveiYiIONWuDcuXw5gxEBIC8+fDlVean/7OJxKi+Ph4unbtSlxcXL79KSkpnDx5Mt/+Bg0aUKtWLZKSkgBISkqiUaNGREVFOc/p3LkzWVlZbN26tcDPO378OFlZWfk2ERGRQBAcDE88AevWmV5Fv/8O//iHGTHKzrY7Os/x+oRo9uzZbNiwgVGjRp11LC0tjZCQECIjI/Ptj4qKIi0tzXnO6clQ3vG8YwUZNWoUERERzi0mJsYNdyIiIuI7GjeGtWtNcuRwmDlFTZqYOUb+yKsTor179zJo0CBmzJhBuXLlSu1zhw4dSmZmpnPbu3dvqX22iIiItyhXzjw+W74catWCX36B9u3h+efB32adeHVClJKSwsGDB7nqqqsoU6YMZcqUYeXKlYwbN44yZcoQFRXFiRMnyMjIyPe+9PR0oqOjAYiOjj6r6izvdd45ZwoNDSU8PDzfJiIiEqiuu85MuO7TB3JzYcQIiI01/Yv8hVcnRB07diQ1NZWNGzc6txYtWtCrVy/nr8uWLcuyZcuc79m+fTt79uwhNjYWgNjYWFJTUzl42jT5JUuWEB4eTsOGDUv9nkRERHxRRIQpzf/0U6hc2XS2vuoqePttkyT5Oodl+VaXgeuuu46mTZvy1ltvAfDwww/z1VdfMXXqVMLDw3n00UcB+PZ/DzlzcnJo2rQp1atX5/XXXyctLY177rmH+++/n5EjRxbpM7OysoiIiCAzM1OjRSIiEvAOHIB+/WDRIvP6hhtgyhSoUcPeuM7kyve3V48QFcXYsWO5+eab6dGjB9dccw3R0dEkJCQ4jwcHB7NgwQKCg4OJjY2ld+/e9OnTh1deecXGqEVERHxXtWpm2Y9334WwMFiyxFSkffyx3ZEVn8+NENlBI0QiIiIF274d7rnHlOkD3H23WSftjAJwWwTUCJGIiIjY57LL4Jtv4IUXTA+jmTPNaNHy5XZH5holRCIiIlIiZcvCyy+bxOiSS+C336BjRxgyBI4dszu6olFCJCIiIm7RujV8/z089JB5PXYstGgBGzfaGlaRKCESERERtylfHiZMgC+/hKgo2LoVWrWC0aMhJ8fu6AqnhEhERETc7qabYMsWuO02OHkSnnnGNHjctcvuyAqmhEhEREQ84sILYc4c06OoYkVYs8askTZlCnhbjbsSIhEREfEYhwPuvRc2bTLroGVnm6aO3bvD77/bHd3/U0IkIiIiHlenDiQmmrlEZcvCvHmmPP/LL+2OzFBCJCIiIqUiOBieeso0cbziCkhPh5tvNlVp2dn2xqaESEREREpVkyawfr3pUwTw3nvQrJlZI80uSohERESk1JUrB2+8AcuWQUwM1KsH0dH2xVPGvo8WERGRQNehA2zeDMePmwnYdlFCJCIiIrbyhoVg9chMREREAp4SIhEREQl4SohEREQk4CkhEhERkYCnhEhEREQCnhIiERERCXhKiERERCTgKSESERGRgKeESERERAKeEiIREREJeEqIREREJOApIRIREZGAp4RIREREAp5Wuy8Cy7IAyMrKsjkSERERKaq87+287/FzUUJUBEeOHAEgJibG5khERETEVUeOHCEiIuKc5zisoqRNAS43N5f9+/dTsWJFHA6HW6+dlZVFTEwMe/fuJTw83K3X9ia6T/8SKPcJgXOvuk//ovs0LMviyJEjVK9enaCgc88S0ghREQQFBVGzZk2PfkZ4eLhf/0+bR/fpXwLlPiFw7lX36V90n5x3ZCiPJlWLiIhIwFNCJCIiIgFPCZHNQkNDefHFFwkNDbU7FI/SffqXQLlPCJx71X36F92n6zSpWkRERAKeRohEREQk4CkhEhERkYCnhEhEREQCnhIiERERCXhKiGw0fvx4ateuTbly5WjdujVr1661OyS3W7VqFd26daN69eo4HA7mzZtnd0geMWrUKFq2bEnFihWpWrUqt956K9u3b7c7LLebMGECjRs3djZBi42NZeHChXaH5XGvvfYaDoeDwYMH2x2K27300ks4HI58W4MGDewOyyP27dtH7969qVKlCmFhYTRq1Ij169fbHZZb1a5d+6z/ng6Hg/j4eLtDc6ucnByGDRtGnTp1CAsLo169egwfPrxIa5YVRgmRTT7++GOGDBnCiy++yIYNG2jSpAmdO3fm4MGDdofmVkePHqVJkyaMHz/e7lA8auXKlcTHx/Pdd9+xZMkSTp48SadOnTh69KjdoblVzZo1ee2110hJSWH9+vV06NCBf/zjH2zdutXu0Dxm3bp1vPfeezRu3NjuUDzmiiuu4MCBA85tzZo1dofkdn/++Sdt27albNmyLFy4kG3btvHGG29QqVIlu0Nzq3Xr1uX7b7lkyRIAbr/9dpsjc6/Ro0czYcIE3nnnHX744QdGjx7N66+/zttvv138i1pii1atWlnx8fHO1zk5OVb16tWtUaNG2RiVZwHW3Llz7Q6jVBw8eNACrJUrV9odisdVqlTJev/99+0OwyOOHDliXXLJJdaSJUusa6+91ho0aJDdIbndiy++aDVp0sTuMDzu6aefttq1a2d3GKVu0KBBVr169azc3Fy7Q3Grrl27Wv369cu3r3v37lavXr2KfU2NENngxIkTpKSkEBcX59wXFBREXFwcSUlJNkYm7pKZmQlA5cqVbY7Ec3Jycpg9ezZHjx4lNjbW7nA8Ij4+nq5du+b7s+qPduzYQfXq1albty69evViz549dofkdvPnz6dFixbcfvvtVK1alWbNmjF58mS7w/KoEydOMH36dPr16+f2hcnt1qZNG5YtW8ZPP/0EwKZNm1izZg1dunQp9jW1uKsN/vjjD3JycoiKisq3Pyoqih9//NGmqMRdcnNzGTx4MG3btuXKK6+0Oxy3S01NJTY2lmPHjlGhQgXmzp1Lw4YN7Q7L7WbPns2GDRtYt26d3aF4VOvWrZk6dSqXXXYZBw4c4OWXX6Z9+/Zs2bKFihUr2h2e2/zyyy9MmDCBIUOG8Oyzz7Ju3Toee+wxQkJC6Nu3r93hecS8efPIyMjg3nvvtTsUt3vmmWfIysqiQYMGBAcHk5OTw4gRI+jVq1exr6mESMTN4uPj2bJli1/OwwC47LLL2LhxI5mZmXz22Wf07duXlStX+lVStHfvXgYNGsSSJUsoV66c3eF41On/om7cuDGtW7fm4osv5pNPPqF///42RuZeubm5tGjRgpEjRwLQrFkztmzZwsSJE/02Ifrggw/o0qUL1atXtzsUt/vkk0+YMWMGM2fO5IorrmDjxo0MHjyY6tWrF/u/pxIiG1x44YUEBweTnp6eb396ejrR0dE2RSXuMHDgQBYsWMCqVauoWbOm3eF4REhICPXr1wegefPmrFu3jv/85z+89957NkfmPikpKRw8eJCrrrrKuS8nJ4dVq1bxzjvvcPz4cYKDg22M0HMiIyO59NJL2blzp92huFW1atXOStovv/xy5syZY1NEnvXrr7+ydOlSEhIS7A7FI5588kmeeeYZ7rrrLgAaNWrEr7/+yqhRo4qdEGkOkQ1CQkJo3rw5y5Ytc+7Lzc1l2bJlfjsXw99ZlsXAgQOZO3cuy5cvp06dOnaHVGpyc3M5fvy43WG4VceOHUlNTWXjxo3OrUWLFvTq1YuNGzf6bTIEkJ2dzc8//0y1atXsDsWt2rZte1YrjJ9++omLL77Ypog8a8qUKVStWpWuXbvaHYpH/PXXXwQF5U9hgoODyc3NLfY1NUJkkyFDhtC3b19atGhBq1ateOuttzh69Cj33Xef3aG5VXZ2dr5/ae7atYuNGzdSuXJlatWqZWNk7hUfH8/MmTP5/PPPqVixImlpaQBEREQQFhZmc3TuM3ToULp06UKtWrU4cuQIM2fOZMWKFSxevNju0NyqYsWKZ83/Kl++PFWqVPG7eWFPPPEE3bp14+KLL2b//v28+OKLBAcH07NnT7tDc6vHH3+cNm3aMHLkSO644w7Wrl3LpEmTmDRpkt2huV1ubi5Tpkyhb9++lCnjn1/z3bp1Y8SIEdSqVYsrrriC77//njfffJN+/foV/6IlrHyTEnj77betWrVqWSEhIVarVq2s7777zu6Q3C4xMdECztr69u1rd2huVdA9AtaUKVPsDs2t+vXrZ1188cVWSEiIddFFF1kdO3a0vv76a7vDKhX+WnZ/5513WtWqVbNCQkKsGjVqWHfeeae1c+dOu8PyiC+++MK68sorrdDQUKtBgwbWpEmT7A7JIxYvXmwB1vbt2+0OxWOysrKsQYMGWbVq1bLKlStn1a1b13ruuees48ePF/uaDssqQVtHERERET+gOUQiIiIS8JQQiYiISMBTQiQiIiIBTwmRiIiIBDwlRCIiIhLwlBCJiIhIwFNCJCIiIgFPCZGIiIgEPCVEIuLXVqxYgcPhICMjw+5QRMSLKSESEVulpaXx6KOPUrduXUJDQ4mJiaFbt275Fj+2g8PhcG7h4eG0bNmSzz//3KVr7N69G4fDwcaNGz0TpIi4jRIiEbHN7t27ad68OcuXL2fMmDGkpqayaNEirr/+euLj4+0OjylTpnDgwAHWr19P27Zt+ec//0lqaqrdYYmIByghEhHbPPLIIzgcDtauXUuPHj249NJLueKKKxgyZAjfffed87w333yTRo0aUb58eWJiYnjkkUfIzs52Hv/111/p1q0blSpVonz58lxxxRV89dVX+T4rJSWFFi1acMEFF9CmTRu2b99+3vgiIyOJjo7m0ksvZfjw4Zw6dYrExETn8UWLFtGuXTsiIyOpUqUKN998Mz///LPzeJ06dQBo1qwZDoeD6667znns/fff5/LLL6dcuXI0aNCAd9991+XfPxFxHyVEImKLw4cPs2jRIuLj4ylfvvxZxyMjI52/DgoKYty4cWzdupVp06axfPlynnrqKefx+Ph4jh8/zqpVq0hNTWX06NFUqFAh3/Wee+453njjDdavX0+ZMmXo169fkWM9deoUH3zwAQAhISHO/UePHmXIkCGsX7+eZcuWERQUxG233UZubi4Aa9euBWDp0qUcOHCAhIQEAGbMmMELL7zAiBEj+OGHHxg5ciTDhg1j2rRpRY5JRNzMEhGxQXJysgVYCQkJLr/3008/tapUqeJ83ahRI+ull14q8NzExEQLsJYuXerc9+WXX1qA9ffffxf6GYBVrlw5q3z58lZQUJAFWLVr17YOHTpU6Ht+//13C7BSU1Mty7KsXbt2WYD1/fff5zuvXr161syZM/PtGz58uBUbG1votUXEszRCJCK2sCyryOcuXbqUjh07UqNGDSpWrMg999zDoUOH+OuvvwB47LHHePXVV2nbti0vvvgimzdvPusajRs3dv66WrVqABw8ePCcnzt27Fg2btzIwoULadiwIe+//z6VK1d2Ht+xYwc9e/akbt26hIeHU7t2bQD27NlT6DWPHj3Kzz//TP/+/alQoYJze/XVV/M9bhOR0qWESERscckll+BwOPjxxx/Ped7u3bu5+eabady4MXPmzCElJYXx48cDcOLECQDuv/9+fvnlF+655x5SU1Np0aIFb7/9dr7rlC1b1vlrh8MB4Hy0VZjo6Gjq169Pp06dmDJlCnfeeWe+JKpbt24cPnyYyZMnk5ycTHJycr64CpI392ny5Mls3LjRuW3ZsiXfvCkRKV1KiETEFpUrV6Zz586MHz+eo0ePnnU8r29QSkoKubm5vPHGG1x99dVceuml7N+//6zzY2JieOihh0hISOBf//oXkydPdmu8rVq1onnz5owYMQKAQ4cOsX37dp5//nk6duzI5Zdfzp9//pnvPXnzjXJycpz7oqKiqF69Or/88gv169fPt+VNwhaR0qeESERsM378eHJycmjVqhVz5sxhx44d/PDDD4wbN47Y2FgA6tevz8mTJ3n77bf55Zdf+Oijj5g4cWK+6wwePJjFixeza9cuNmzYQGJiIpdffrnb4x08eDDvvfce+/bto1KlSlSpUoVJkyaxc+dOli9fzpAhQ/KdX7VqVcLCwli0aBHp6elkZmYC8PLLLzNq1CjGjRvHTz/9RGpqKlOmTOHNN990e8wiUkR2T2ISkcC2f/9+Kz4+3rr44outkJAQq0aNGtYtt9xiJSYmOs958803rWrVqllhYWFW586drQ8//NACrD///NOyLMsaOHCgVa9ePSs0NNS66KKLrHvuucf6448/LMv6/0nVeedalmV9//33FmDt2rWr0LgAa+7cufn25ebmWg0aNLAefvhhy7Isa8mSJdbll19uhYaGWo0bN7ZWrFhx1vsmT55sxcTEWEFBQda1117r3D9jxgyradOmVkhIiFWpUiXrmmuuKdYEcxFxD4dluTCzUURERMQP6ZGZiIiIBDwlRCIiIhLwlBCJiIhIwFNCJCIiIgFPCZGIiIgEPCVEIiIiEvCUEImIiEjAU0IkIiIiAU8JkYiIiAQ8JUQiIiIS8JQQiYiISMD7P/PEFs+Y+DByAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(cash_rate, house_prices, 'ro')\n",
    "plt.xlabel(\"Cash Rate\")\n",
    "plt.ylabel(\"House Prices\")\n",
    "\n",
    "# Here we create x values for \"each hour of study\" so we can plot the results\n",
    "x_model = np.arange(0, 8, 0.25)\n",
    "y_model = x_model * beta[1] + beta[0]  # Use our learned model to fit the line of best fit.\n",
    "\n",
    "plt.plot(x_model, y_model, 'b-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predicting the actual scores, we can see our trained model is better than the intuited model from earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The squared error is 438881.87\n"
     ]
    }
   ],
   "source": [
    "predicted_prices = cash_rate * beta[1] + beta[0]\n",
    "squared_error = np.sum((predicted_prices - house_prices) ** 2)\n",
    "print(\"The squared error is {:.2f}\".format(squared_error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercises\n",
    "\n",
    "1. Remove the values with a cash rate above 6.0 and rerun the analysis. What $\\beta$ values do you obtain?\n",
    "2. Implement and run the following algorithm:\n",
    "\n",
    "\n",
    "```\n",
    "loop 100 times:\n",
    "    obtain a sample of 50% of the data points\n",
    "    compute the beta values\n",
    "    \n",
    "average all beta values obtained during the loop as the final beta values\n",
    "```\n",
    "\n",
    "This is an *ensemble* learner - it is great for removing the effect of outliers, as seen in the above data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The squared error is 219662.82\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVoUlEQVR4nO3de3zO9f/H8ce1sZnDNhabw3JIhVKGMKeISKKiQkgpOkwZ31Q6kHIWlUOJ+jqUQ8qh0g8NQ6JhTKh0UoRRDptDhu3z++OdfY1tdm3Xtc+1Xc/77Xbd5vpc730+r8+3b11Pn/fJYVmWhYiIiIgX87G7ABERERG7KRCJiIiI11MgEhEREa+nQCQiIiJeT4FIREREvJ4CkYiIiHg9BSIRERHxekXsLqAgSEtL48CBA5QqVQqHw2F3OSIiIpIDlmVx4sQJKlSogI9P9s+AFIhy4MCBA4SHh9tdhoiIiOTCvn37qFSpUrZtFIhyoFSpUoD5HzQwMNDmakRERCQnkpOTCQ8PT/8ez44CUQ5c6CYLDAxUIBIRESlgcjLcRYOqRURExOspEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8ngKRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5PgchmGzbA33/bXYWIiIh3UyCy0b590KED1KkDX39tdzUiIiLeS4HIRidPQtmysH8/tGwJI0dCWprdVYmIiHgfBSIb1awJW7ZAz56QmgovvQR33AGHD9tdmYiIiHdRILJZyZIwaxb8978QEAAxMXDzzRAba3dlIiIi3kOByAM4HPDII7B5M9SqBYmJ0Lo1DBtmnhyJiIiIeykQeZAbbjChqHdvM5bo1VehTRs4eNDuykRERAo3BSIPU7w4fPABzJ4NJUrA6tVmFlpMjN2ViYiIFF4KRB6qZ08z4Lp2bTPIum1bePllOH/e7spEREQKHwUiD1ajBsTFweOPg2XBiBFw223w5592VyYiIlK4KBB5uIAAmDoV5s2DUqXMAo516sCyZXZXJiIiUngoEBUQXbtCfDxERMCRI3DnnfD883DunN2ViYiIFHwKRAXItdeavc+iosz7sWPh1lth71576xIRESnoFIgKmGLFYPJk+PRTCAqCjRtNF9rnn9tdmYiISMGlQFRAde4MW7fCLbfAsWNw990wcCCcPWt3ZSIiIgWPAlEBVq0arF8P0dHm/ZtvQrNmsGePrWWJiIgUOApEBVFqKqxZA/Pm4bdhDW++kcpnn0Hp0rBpkxl4vWiR3UWKiIgUHApEBc2iRVClCrRsCQ8+aH5WqULH84vYtg0aNYKkJNOl9vTTkJJid8EiIiKeT4GoIFm0CO677/KVGffvh/vuo3L8Itatg+eeM4cnT4bGjeGXX/K/VBERkYJEgaigSE2F/v3NktWXunAsOpqiPqmMGQNffgkhIWbgdd26sGCBi2r4t6uONWvMexERkUJAgaig+Prr7PfssCzYt8+0wyzcmJAATZvCiRPQpQs88QT8808ur59FV50GK4mISGGgQFRQHDzodLtKlSA2Fl56CRwOeO89M8Zo924nr32FrjqFIhERKegUiAqK8uVz1a5IERg+HFasgHLl4LvvoF49mDMnh9fNYVddtt1n6moTEREPp0BUUDRrZh75OByZf+5wQHi4aZeJ2283XWgtW8KpU9CjBzz2GJw+fYXrOtlVdxl1tYmISAFgayBat24dHTp0oEKFCjgcDpYsWZLhc8uyGDJkCOXLlycgIIDWrVvz888/Z2hz9OhRunfvTmBgIMHBwTz66KOcPHkyQ5vvvvuOZs2aUaxYMcLDwxk7dqy7b831fH3h7bfNny8NRRfev/WWaZeF8uUhJgZefdX8ygcfQIMG8P332Vw3F1116dTVJiIiBYStgejUqVPcfPPNTJkyJdPPx44dy8SJE5k6dSpxcXGUKFGCtm3bcubMmfQ23bt3Z9euXcTExLB06VLWrVtH37590z9PTk6mTZs2VK5cmfj4eMaNG8err77KtGnT3H5/Ltepk9nErGLFjMcrVTLHO3W64il8fWHoUFi1CsLCYNcus/3HzJlZ/EIuu+pc0tV24TzqbhMREXezPARgLV68OP19WlqaFRYWZo0bNy792PHjxy1/f39r3rx5lmVZ1vfff28B1ubNm9PbLFu2zHI4HNb+/fsty7Ksd955xypdurSVkpKS3ub555+3rr/++hzXlpSUZAFWUlJSbm/Ptc6ft6zYWMuaO9f8PH8+V6dJTLSs22+3LJNQLOuhhyzrxIlMrlWpkmU5HP9rePHL4bCs8PDLa4iNzbz9pa/Y2KwLXLjQXPvi9pUqmeMiIiJX4Mz3t8eOIdqzZw+JiYm0bt06/VhQUBANGzZk48aNAGzcuJHg4GDq16+f3qZ169b4+PgQFxeX3qZ58+b4+fmlt2nbti27d+/m2LFjmV47JSWF5OTkDC+P4usLLVpAt27mZzbdZNkJDYXly2HECPDxgdmzzdOiHTsuuVZuuury0tUGee9u05MlERFxgscGosTERABCQ0MzHA8NDU3/LDExkXLlymX4vEiRIpQpUyZDm8zOcfE1LjVq1CiCgoLSX+Hh4Xm/IQ/l4wMvvmgyQ8WK8OOPZlzR9OkX9Xblpqsut11tkPfuNg3kFhERJ3lsILLT4MGDSUpKSn/t27fP7pLcrlkzMwutXTs4cwb69jVZIv3hWKdO8PvvZmGjuXPNzz17sh63lJdZcXmZ2aaB3CIikgseG4jCwsIAOHToUIbjhw4dSv8sLCyMw4cPZ/j8/PnzHD16NEObzM5x8TUu5e/vT2BgYIaXN7jqKli6FMaONT1g8+ebNYu2bfu3gTNddXmZFZfb7jZXDeQWERGv47GBqGrVqoSFhbFq1ar0Y8nJycTFxREZGQlAZGQkx48fJz4+Pr3N6tWrSUtLo2HDhult1q1bx7lz59LbxMTEcP3111O6dOl8upuCw8cHBg0yD1/Cw83GsI0awTvvZJ4zspXbWXG57W7L65pJIiLitWwNRCdPniQhIYGEhATADKROSEhg7969OBwOoqOjGT58OJ9//jk7duzgoYceokKFCtxzzz0A1KxZkzvuuIM+ffqwadMmvvnmG/r160fXrl2pUKECAA8++CB+fn48+uij7Nq1i48//pi3336bgQMH2nTXBUNkpOlC69gRzp6FqCh44AE4ftzJEznb1Qa5727L60BuERHxXvkw6y1LsbGxFnDZq1evXpZlman3r7zyihUaGmr5+/tbrVq1snbv3p3hHEeOHLG6detmlSxZ0goMDLQeeeQR68Qlc8e3b99uNW3a1PL397cqVqxojR492qk6PW7afT5KS7OsCRMsq2hRM+u9alXL2rQpHy68cKGZ0n/pdP8LxzKber9yZc6m+q9cmQ83ICIidnPm+9thWU53hHid5ORkgoKCSEpK8prxRJfatAm6dDEPe4oWNeOM+vfP+iGOSyxaZC5ycTdYeLgZe5TZE6ZVq+CiZRqytHIltGrlsjJFRMQzOfP97bFjiMSzNGhgBld36gTnzsGAAXDvvXD0qBsv6mx3WxbLKOS6nYiIeA0FIsmx4GAzFnrSJPDzg88+g4gI+PZbN17UmZltf/2Vs3PmtJ2IiHgNBSJxisMB/frBxo1wzTWwd68Z2zxuHKSl2VxcSIhr24mIiNdQIJJcqVsXtm4144rOn4fnnoMOHeDvv20s6sgR17YTERGvoUAkuRYYaLYKmzoV/P3h//4P6tSxcZmfsmVd205ERLyGApHkicMBjz8OcXFw3XVmh4yWLWHkSBu60LJYeTzX7URExGsoEIlL3HwzxMdDjx5mZ4yXXjL7ol2ys4qIiIhHUiASlylZEmbPhg8+gIAA+Oor04W2Zk0+FZDT9KWUJiIil1AgEpdyOKB3b9i8GWrWNLtktGoFw4blw56q5cq5tp2IiHgNBSJxixtuMKHokUfMWKJXX4U2bbQmooiIeCYFInGbEiXgv/813WglSsDq1Was0cqVbrqgusxERCSXFIjE7Xr2hC1boHZtk0XatIGXXzbrF7lU+fKubSciIl5DgUjyRY0aZmp+375my/kRI8zYov37XXiRxo2z39oDzOeNG7vwoiIiUhgoEEm+CQiA994zizmWLAnr1plZaMuXu+gCGzZceeR2aqppJyIichEFIsl3XbuabT8iIsxWH+3awQsvwLlzeTzxwYOubSciIl5DgUhsce215kFNVJR5P2aM2cx+7948nFRjiEREJJcUiMQ2xYrB5Mnw6acQFGQCUkQEfPFFLk+oMUQiIpJLCkRiu86dTRfaLbfA0aPQsSP85z9w9qyTJ9IYIhERySUFIvEI1arB+vUQHW3eT5gAzZrB7787cRKNIRIRkVxSIBKP4ecHb74Jn30GpUvDpk2mC23x4hyeQGOIREQklxSIxON07AjbtkGjRnD8OHTqBM88AykpV/jFZs2gUqXs24SHm3YiIiIXUSASj1S5slmn6LnnzPtJk6BJE/j112x+ydcXunXL/sRdu1554LWIiHgdBSLxWEWLmun4X34JISEQH2+60BYsyOIXUlPNqo/ZmT//ygOvRUTE6ygQice7805ISICmTeHECejSBZ58Es6cuaTh11/Dn39mf7J9+0w7ERGRiygQSYFQqRLExsKLL4LDAVOnmjFGP/10USPNMhMRkVxSIJICo0gRsyns8uVQtixs3w5168KcOf820CwzERHJJQUiKXDatDFhqEULOHUKevSAxx6D0/U0y0xERHJHgUgKpPLlYeVKGDrUdKF98AE0bOzLD236Z/+LmmUmIiKZUCCSAsvXF1591QSjsDDYuRPqz3iKWTyU9S9plpmIiGRCgUgKvNtuM7PQWtc7ymmrOA8zi17M5BTFL2+sWWYiIpIJBSIpFEJDYcWAFQznJXxIZTa9qM8WdnDj5Y01y0xERC6hQCSFhk/F8rzESGJpSQX28yM1acAmpvMY1sUNNctMREQuoUAkhce/e5k152sSqMMdLOMMAfRlOt2ZwwlKapaZiIhkSoFICo+L9jIry998SXvG8By+nGceD1KPeBJaDtAsMxERuYwCkRQel+xl5oPFc4xjHc0JZy8/cx2NZj/JO5PTsKxsziMiIl5HgUgKjyz2MmvMRrYRQQc+J4ViRD3twwMPQFKSDTWKiIhHUiCSwiOb2WMhHOUz7mYCAyjim8ann5ptP7Zsycf6RETEYykQSeFxhdljDmAAb/HNpG1UqQK//QaNG8Pbb5OxCy01FdasMd1va9ZoIUcRES+gQCSFx7+zzHA4Mv/c4YDwcBr0rcO2bXDvvXDuHERHQ6dOcOwYsGgRVKkCLVvCgw+an1WqmOMiIlJoKRBJ4eHrax73wOWh6ML7t94CX1+Cg2HhQpg0Cfz8YMkSiLj+FN92Hnf5OKT9++G++xSKREQKMQUiKVw6dYJPP4WKFTMer1TJHO/UKf2QwwH9+sGGDXDNNRZ//FWCZqzjDf5DGhcFqgv9adHR6j4TESmkHJalCchXkpycTFBQEElJSQQGBtpdjuREaqqZdXbwoBlb1KxZtusPJS39mr4dDrCALgC0ZykzeZirOJKxYWwstGjhxsJFRMRVnPn+LpJPNYnkL19fp4JL0Ik/mc+DtCSWaN7iS+4igm3MoxtN+eZ/DbUPmohIoaQuMxGA8uVxAE/wHnE05Dp28yfhtGANo3jhf11o2gdNRKRQUiASgQwz1G7mO7ZQn+58RCpFeJFRtGMZhyvUyX4fNE3XFxEpsBSIROCyGWqlOMmH9OQDehPAab6iLXX+2cCar7MYh6Tp+iIiBZoCkcgFl8xQcwC9mcGm0I7UrJTMwWMBtGoFr712ycOfRYvMtHxN1xcRKbA0yywHNMvMy2QyQ+3UGV/69YOZM02T226DOXMgrGyqeRKUyR5qgJnbX6kS7NmT7Sw3ERFxPWe+v/WESORSF2aodetmfvr6UqIEzJgBs2ZB8eKwejXUqQOr3tqRdRgCs4bRvn0mYImIiMdSIBJxwkMPQXw83HgjHDoEtw+6mSEM4zxXePqj6foiIh5NgUjESTVqwKZN0KcPWJaD1xlCK1ZxgGym5Gu6voiIR1MgEsmFgACYNg3mfpRGScdJ1nErN7Od5bS9vHF4ePbT9UVExHYKRCJ50K27D/G936EO2/ibsrRjOYMZybmLF4Hv2lUDqkVEPJwCkUhepKZy3YpJbCSSp5gCwGgG05JY9lHJtJk/X4s0ioh4OAUikbz4+mv480+KkcIU+rGA+wkkiW9oSh0SWEp7zTITESkAFIhE8uKS2WP38ynbiKA+mzlKCB1YyrOM4+y+QzYVKCIiOaFAJJIXmcweq8Ye1tOU/rwFwHiepfmodvz+e/6WJiIiOadAJJIXF20KezF/zvIWA1jMvQQ7koj7IZCICFiyxJ4yRUQkewpEInlxYVPYLHbAuYclJLyzgYYN4fhxuPde6N8fUlLyt0wREcmeApGIm1Uu9w9ffw3PPmveT5wITZrAr7/aW5eIiPyPApFIXqSmmkc+WXE4IDqaoj6pjBsHS5dCSIjZ/qNuXfjkk/wrVUREsubRgSg1NZVXXnmFqlWrEhAQwDXXXMPrr7+OdVH3hGVZDBkyhPLlyxMQEEDr1q35+eefM5zn6NGjdO/encDAQIKDg3n00Uc5efJkft+OFEb/TrvP0iWbu7ZvDwkJ5glRcjI88AA89RScOZM/5YqISOY8OhCNGTOGd999l8mTJ/PDDz8wZswYxo4dy6RJk9LbjB07lokTJzJ16lTi4uIoUaIEbdu25cxF3zDdu3dn165dxMTEsHTpUtatW0ffvn3tuCUpbHK6aetF7SpVgjVrYPBg8/7dd6FRI/jpJ9eXJyIiOeOwrCxGg3qAu+66i9DQUD744IP0Y507dyYgIICPPvoIy7KoUKEC//nPf3j23wEaSUlJhIaGMnPmTLp27coPP/xArVq12Lx5M/Xr1wdg+fLl3Hnnnfz5559UqFDhsuumpKSQctGo1+TkZMLDw0lKSiIwMNDNdy0Fypo10LLlldvFxkKLFpcdXrECevaEv/6CkiXhvffgwQddXqWIiFdKTk4mKCgoR9/fHv2EqHHjxqxatYqf/v2r8/bt21m/fj3t2rUDYM+ePSQmJtK6dev03wkKCqJhw4Zs3LgRgI0bNxIcHJwehgBat26Nj48PcXFxmV531KhRBAUFpb/Cw8PddYtS0GUx7T6dw5Ht5q5t25outBYt4ORJ6N4d+vSB06fdVrGIiGTCowPRCy+8QNeuXalRowZFixYlIiKC6OhounfvDkBiYiIAoaGhGX4vNDQ0/bPExETKlSuX4fMiRYpQpkyZ9DaXGjx4MElJSemvffv2ufrWpLC4MO0eLg9FF96/9Va2m7tWqAArV8KQIeZX3n8fGjaEH35wT8kiInI5jw5ECxYsYM6cOcydO5etW7cya9Ys3njjDWbNmuXW6/r7+xMYGJjhJZKlTp3MnHqfS/518vExxzt1uuIpfH1h2DCIiYHQUNi5E+rXBzf/X11ERP7l0YFo0KBB6U+JateuTc+ePRkwYACjRo0CICwsDIBDhzLuE3Xo0KH0z8LCwjh8+HCGz8+fP8/Ro0fT24jkyaJF8MYbl+9on5pqji9alONTtWplutBatTLdZg8/bF6nTrmyYBERuZRHB6LTp0/jc8nfun19fUlLSwOgatWqhIWFsWrVqvTPk5OTiYuLIzIyEoDIyEiOHz9OfHx8epvVq1eTlpZGw4YN8+EupFC7sA5RdnMToqMvD0vZCAszg61ff908ZJo1yzwt2rkz7+WKiEjmPDoQdejQgREjRvDll1/y+++/s3jxYiZMmMC9994LgMPhIDo6muHDh/P555+zY8cOHnroISpUqMA999wDQM2aNbnjjjvo06cPmzZt4ptvvqFfv3507do10xlmIk5xch2inPL1hZdfhtWrzRijH3+EW24x44s8d16oiEjBVcTuArIzadIkXnnlFZ566ikOHz5MhQoVePzxxxkyZEh6m+eee45Tp07Rt29fjh8/TtOmTVm+fDnFihVLbzNnzhz69etHq1at8PHxoXPnzkycONGOW5LCJhfrEDnj1ltNF9pDD8Hy5WYGWmwsTJ0KpUrl6pQiIpIJj16HyFM4s46BeJk8rkOUU2lpMG4cvPSS6X277jr4+GOoUyfXpxQRKfQKzTpEIh4vj+sQpUtNNeFq3jzz85IxRz4+8PzzsHatudxPP5nVrd99V11oIiKuoEAkkhcuWIeIRYugShXzpOnBB83PKlUynZ3WpInpQrvrLkhJMfugdekCSUkuuBcRES+mQCSSV506waefQsWKGY9XqmSOZ7cO0aJFcN99lw/M3r/fHM8kFIWEwOefw/jxUKQIfPIJ1K0LF02kFBERJ2kMUQ5oDJHkSGqqmU128CCUL2+6ybJ7MpSaap4EZTVLzeEwoWrPnizPExdnnhD98QcULWqWPXr66ax78EREvInGEInYwdfXDJzu1s38zC4MgUum7DdsCNu2wb33wrlzZkmkzp3h2LFc3YGIiNdSIBKxi4um7JcuDQsXwsSJ4OcHixdDRIR5eiQiIjmjQCRil/LlXdbO4TBdZRs2QLVqpgutaVMzzkid4iIiV6ZAJGIXV0zZv2S6fr06qWzdCvffD+fPm71lO3aEI0fccgciIoWGApGIXfI6ZT+L6fpBqxbx8cdmjSJ/f1i61HShffONu25ERKTgUyASsVNup+xfYbq+Y/EinngCvv0Wrr3WjM2+9VYYPdqsei0iIhlp2n0OaNq9uJ0zU/adnK5/4gQ88QTMnWs+vuMOmD0bypZ1y52IiHgMTbsXKWicmbLv5HT9UqXgo4/g/fehWDGzSWydOmYbEBERMRSIRAqaXEzXdzjg0Udh82aoUQMOHIDbboPXX79s2zQREa+kQCRS0ORhuv6NN8KWLdCrlxlLNGQItG0LiYkurlFEpIBRIBIpaPI4Xb9ECZg507yKF4dVq0wX2qpV7ipYRMTzKRCJFDQXputnNR/CsrKfrv+vXr1MF9qNN8KhQ3D77TB0qLrQRMQ7KRCJeLFatcwWH489ZnLUa69Bq1ZmjJGIiDdRIBIpaFJTzS6uWXE4IDo6x496iheH6dNhzhwoWdLMPqtTB1ascEm1IiIFggKRSEHj5LT7nHrwQYiPh5tvhr/+MusVDR5stgARESnsFIhECppcTLvPqeuuM6tbP/mkeT96tFkWad8+p08lIlKgKBCJFDR5mHafE8WKwTvvwMcfQ2Cg2QOtTh348stcnU5EpEBwOhDNmjWLLy/6L+Nzzz1HcHAwjRs35o8//nBpcSKSiQvT7rOTzbT7nHrgAdi6FerVg6NH4a674Nln4dy5PJ1WRMQjOR2IRo4cSUBAAAAbN25kypQpjB07lquuuooBAwa4vEARuYSvr9niIztdu15x2n1OXHONeUL0zDPm/fjxJmf9/nueTy0i4lGcDkT79u2jevXqACxZsoTOnTvTt29fRo0axddODuIUkVxITYV587JvM3++yxYU8vc3yx4tXgzBwWaafkQELFniktOLiHgEpwNRyZIlOXLkCABfffUVt99+OwDFihXjn3/+cW11InK5K80yg1zNMruSe+6BbdugYUM4fhzuvdfM/k9JcellRERs4XQguv3223nsscd47LHH+Omnn7jzzjsB2LVrF1WqVHF1fSJyKTfOMruSKlVMznr2WfN+4kRo0gR++83llxIRyVdOB6IpU6YQGRnJX3/9xcKFCwkJCQEgPj6eblca1yAieefmWWZXUrQojBsHX3wBZcqYtYsiIuDTT91yORGRfOGwrKw2RJILkpOTCQoKIikpicDAQLvLEW+Xmmoe1WTXbRYeDnv2uGRgdXb27TPju7/5xrx/8kmYMMFM3RcRsZsz39+5Wofo66+/pkePHjRu3Jj9+/cD8OGHH7J+/frcnE5EnJGPs8yuJDwcYmPNitYA774LjRrBTz+5/dIiIi7ldCBauHAhbdu2JSAggK1bt5Ly74jKpKQkRo4c6fICReQS+TzL7EqKFoWRI2H5cihbFrZvN2sXzZ2bL5cXEXEJpwPR8OHDmTp1KtOnT6do0aLpx5s0acLWrVtdWpyIZMKmWWZX0rYtJCTArbfCyZPQvTv06QOnT+drGSIiueJ0INq9ezfNmze/7HhQUBDHjx93RU0ikh0bZ5ldSYUKsHIlDBkCDge8/76Zpv/DD/leioiIU5wORGFhYfzyyy+XHV+/fj3VqlVzSVEikg2bZ5ldSZEiMGwYxMRAaCjs3An168OsWbaUIyKSI04Hoj59+tC/f3/i4uJwOBwcOHCAOXPm8Oyzz/LkhS2yRcR9mjWDf5e7yFJISJ73MsurVq1MF1qrVqbb7OGHzevUKVvLEhHJVBFnf+GFF14gLS2NVq1acfr0aZo3b46/vz/PPvssTz/9tDtqFJECKiwMVqyAUaNg6FDzlGjTJliwAG680e7qRET+J9frEJ09e5ZffvmFkydPUqtWLUqWLOnq2jyG1iESj7JmDbRseeV2sbHQooW7q8mxtWvhwQfhwAGzTtHkydC7txlrJCLiDm5dhygpKYmjR4/i5+dHrVq1aNCgASVLluTo0aMkJyfnumgRySEPHlSdnVtvNV1od9wBZ87AY49Bz55w4oTdlYmI5CIQde3alfnz5192fMGCBXTt2tUlRYlINjx8UHV2ypaFL7+E0aPNupFz5pgB19u3212ZiHg7pwNRXFwcLTN5XN+iRQvi4uJcUpSIZKNx4yuvQu3ra9p5IB8feP5504VWqZJZ1bphQ5g6FbSRkIjYxelAlJKSwvnz5y87fu7cOf755x+XFCUi2diw4cqrUKemmnYerEkT04V2112QkmL2QevaFZKS7K5MRLyR04GoQYMGTJs27bLjU6dOpV69ei4pSkSy4aoxRKmpZoD2vHnmZz5t9XGxkBD4/HMYP96sX7Rggdn2Iz4+30sRES/n9LT74cOH07p1a7Zv306rVq0AWLVqFZs3b+arr75yeYEicglXjCFatAj698+4BUilSvD229CpU97qc5LDAQMHmidGXbrAr7+a3r433oB+/TQLTUTyh9NPiJo0acLGjRsJDw9nwYIFfPHFF1SvXp3vvvuOZjYvBCfiFZo1M+Elq6TgcJht6LP693HRIrjvvsv3Q9u/3xxftMi19eZQw4awbRvccw+cPQvPPAOdO8OxY7aUIyJeJtfrEHkTrUMkHudCqIGMI5EvhKRPP838SU9qKlSpkvXmsA6HCVt79mQ9cDs11Wwce/CgeQrVrNmVB3k7wbJg0iR49lk4d86U+/HH0KCByy4hIl7C5esQXby+UHJycrYvEckHnTqZ0FOxYsbjlSplHYbABJmswhCYNLJvn2mXmUWLTEJp2dKsstiypXnvwqdKDod5OrRhA1SrBr//brrTJkzQLDQRcZ8cPSHy9fXl4MGDlCtXDh8fHxyZPKq3LAuHw0GqDQMz3U1PiMRjOfu0Zt48E2SuZO5c6NYt47ELT6Uu/U/GlZ5K5UFSEvTpA598Yt536AAzZ0KZMi69jIgUUs58f+doUPXq1asp8+9/gWJjY/NeoYi4hq+vc9tz5HZAdmqqGYSd2d+fLMuEouhouPtul3afBQWZ7rKWLWHAAPjiC6hTB+bP99hllkSkgHJqDNH58+cZOXIkvXv3plKlSu6sy6PoCZEUGhfGEO3fn3m4yWoMkQfsn5aQAA88AD//bEobMQIGDTILPYqIZMZte5kVKVKEcePGZbowo4gUAL6+Zmo9XD5L7cL7t966/CmPB+yfVqeOWZ+oWzeT6154Adq3h7/+ctslRcSLOP13q9tuu421a9e6oxYRyQ+5GZBdrlzOzp3TdrlUqpTZ/2z6dChWDJYvN0Fp3Tq3XlZEvIDTCzO2a9eOF154gR07dlCvXj1KlCiR4fOOHTu6rDgRcZNOncx4HzdOn3e5fweQOw4e5LHq5Wm4sRkPdPPlxx9Nb96wYTB4sGffgoh4LqfXIfLJpsNes8xECqm8zE5zhSxW1j45ejJPrbibDz80h1q3ho8+gtBQ15cgIgWP28YQAaSlpWX5KoxhSERwzXYhuZXNytole97L7HsWMWMGFC8OK1fCzTfD6tWuL0NECjennhD9/vvvxMTEcO7cOW699VZuuOEGd9bmMfSESLxebmenueq62S0mGR4Oe/bw/W5fHngAdu0y5bzyCgwZoi40EW/mlidEsbGx3HDDDTz++OP069ePiIgIPvroozwXKyIFQG5np+XVlVbWhvSVtWvVgk2b4NFHTWZ77TXThXbggGtLEpHCKceB6JVXXuH2229n//79HDlyhD59+vDcc8+5szYR8SS53S4kL/bvd6pd8eLw/vtmHFGJEmb5pDp14KuvXF9atlJTzcXnzTM/NZxAxOPlOBDt3LmTkSNHUr58eUqXLs24ceM4fPgwR44ccWd9IuJJOnUym4vFxpoB1LGxppvMHWEIcr7I0CXtuneHrVvNeKK//oK2beHFFyFfllDLh/3eRMT1chyIkpOTueqqq9LfFy9enICAAJKSktxSmIh4qAvbhXTrZn66c5BO2bK5bnfddbBxIzzxhHk/apTJJlfqgcuTbAaAc999CkUiHsypdYhWrFhBUFBQ+vu0tDRWrVrFzp07049pHSIRcZlLu+ecbBcQAO++a4LQY4/B+vWmC23WLLPKtUtdab83cMt+byLiGjmeZZbd+kPpJ9M6RCLiSk7MMrtSyPjlF+jSxXSlATz7LIwcCUWLuqhWu/d7+3fhygKz0KZIPnDLLLPs1h9y5zpE+/fvp0ePHoSEhBAQEEDt2rXZsmVL+ueWZTFkyBDKly9PQEAArVu35ueff85wjqNHj9K9e3cCAwMJDg7m0Ucf5eTJky6vVURc7MLstktntl3gcOR4dlv16rBhAzz9tHn/xhvQvDn88YeLanVyALhLadySSJ559D7Rx44do0mTJhQtWpRly5bx/fffM378eEqXLp3eZuzYsUycOJGpU6cSFxdHiRIlaNu2LWfOnElv0717d3bt2kVMTAxLly5l3bp19O3b145bEhFnXZjdVqlSxuPh4U7PbvP3h4kTYeFCCAqCb781XWhLlrigzlwOAM8zjVsScQ3Lgz3//PNW06ZNs/w8LS3NCgsLs8aNG5d+7Pjx45a/v781b948y7Is6/vvv7cAa/Pmzeltli1bZjkcDmv//v05qiMpKckCrKSkpFzeiYjk2fnzlhUba1lz55qf58/n6XS//WZZDRpYlhngY1n9+1tWSkoeTvjRR/87WXavjz7KU90ZnD9vWZUqZX+98PA8/28lUlA58/3t0U+IPv/8c+rXr8/9999PuXLliIiIYPr06emf79mzh8TERFq3bp1+LCgoiIYNG7Jx40YANm7cSHBwMPXr109v07p1a3x8fIiLi8v0uikpKSQnJ2d4iYjNXDy7rWpVM+Rm4EDz/u23oUkT+O23XJ4wjwPAc8WJhStFJHseHYh+++033n33Xa699lpWrFjBk08+yTPPPMOsWbMASExMBCD0kp0cQ0ND0z9LTEykXLlyGT4vUqQIZcqUSW9zqVGjRhEUFJT+Cg8Pd/WtiYgH8POD8ePh88+hdGnYsgUiIkxPnNOaNbu8W+9S4eGmnavYOW5JpJDx6ECUlpZG3bp1GTlyJBEREfTt25c+ffowdepUt1538ODBJCUlpb/27dvn1uuJiL06dICEBGjcGJKT4f77ISoKLhqKeGUuHACeY3aNWxIphHIViI4fP87777/P4MGDOXr0KABbt25lv4v/FlK+fHlq1aqV4VjNmjXZu3cvAGFhYQAcOnQoQ5tDhw6lfxYWFsbhw4czfH7+/HmOHj2a3uZS/v7+BAYGZniJSOF29dVm5vzzz5v377wDkZFwyaTV7LlwAHiOhIS4tp2IF3M6EH333Xdcd911jBkzhjfeeIPjx48DsGjRIgYPHuzS4po0acLu3bszHPvpp5+oXLkyAFWrViUsLIxVq1alf56cnExcXByRkZEAREZGcvz4ceLj49PbrF69mrS0NBo2bOjSekWkYCtaFEaPhv/7P7jqKvPUqG5dsyVZjuV1exNn9kHL6dZJ2mJJ5MqcHbHdqlUra9CgQZZlWVbJkiWtX3/91bIsy/rmm2+sypUrO3u6bG3atMkqUqSINWLECOvnn3+25syZYxUvXtz66KJZGqNHj7aCg4Otzz77zPruu++su+++26patar1zz//pLe54447rIiICCsuLs5av369de2111rdunXLcR2aZSbiff7807KaNfvfZK0+fSzr9Gk3X3ThwstnjVWqZI5nZvbsnM1smz3bzYWLeCZnvr+dDkSBgYHWL7/8YllWxkD0+++/W/7+/s6e7oq++OIL68Ybb7T8/f2tGjVqWNOmTcvweVpamvXKK69YoaGhlr+/v9WqVStr9+7dGdocOXLE6tatm1WyZEkrMDDQeuSRR6wTJ07kuAYFIhHvdO6cZb38smU5HCZX1K5tWT/84KaLLVz4vwtd/HI4zCuzUPTmmzkLRG++6aaiRTybM9/fOd6644Jy5cqxYsUKIiIiKFWqFNu3b6datWrExMTQu3fvQjkAWVt3iHi3mBjo0QMOH4bixc3+aA895MILXGmLEofDjEu6dIuSOXNMYVfy0UfQvbtLShUpSNyydccFHTt25LXXXuPcuXOA2b9s7969PP/883Tu3Dl3FYuIeKJ/x/Pc/vc8Et7ZQMsWFqdPQ69e8MgjcOqUi65zpfWELCvz9YTsWPtIpJByOhCNHz+ekydPUq5cOf755x9uvfVWqlevTqlSpRgxYoQ7ahQRyX+X7A9W/r4mxPxchWFdvsfHB2bOhAYNYNcuF1zr4MHctbNj7SORQqqIs78QFBRETEwM33zzDdu3b+fkyZPUrVs3w2rRIiIeIze7wF/YH+ySEQW+B/YxZMGNNH91DQ9Obc7338Mtt8CkSdC7d9ZLEF3RJYvH5rjdhbWPMqkVcM/aRyKFlNNjiDJz/PhxgoODXVCOZ9IYIpECatEi6N8/Y3dUpUomRGQ1DT6H43kOx+2h58O+fPWVOdy9uxlbVKpULupctQpy8pfKlSuhVavLj2d2n+HhJgy5eu0jkQLErWOIxowZw8cff5z+/oEHHiAkJISKFSuyfft256sVkYLFmXVy7JTbXeBzOJ6n3O6vWbYMRo40D2DmzIH69SFX/xm8ZPFYp9vlde0jEXE+EE2dOjV9b6+YmBhiYmJYtmwZ7dq1Y9CgQS4vUEQ8yCXjamjZ0rzPKlxcKr/CVGqqeWKS2QPwC8eiozO/vhPjeXx8YPBgcysVK8JPP0HDhvDeu2lYsWtyfp/ly+fsmtm1c/HmtyLexulAlJiYmB6Ili5dygMPPECbNm147rnn2Lx5s8sLFBEPsWgRdO58+dOTP/80x68UivIappyR21lbkKtw0rSpWdX6zjshJQWeeMqHrrcdIvnBx3N2nxcGR2e3D5oGR4u4ldOBqHTp0ulrDS1fvjx9MLVlWaR66qNzEcmb1FTo2zf7Nn37Zv0kJLfdV7mV21lbkOtwctVV8MUjixjHIIpwjgV0oS5biafule/zwuDoC+e/9HqgwdEibuZ0IOrUqRMPPvggt99+O0eOHKFdu3YAbNu2jerVq7u8QBHxAGvWXHk/rCNHTLtL5aX7Krfy0gWV23CSmorPgP48yxt8TTOu5g9+pTqN2cBk6ylzq9nd54WNYS9dM6hSJfdsDCsiGTgdiN5880369etHrVq1iImJoWTJkgAcPHiQp556yuUFiogHyCzo5LRdXrqvciuvXVC5CScX3Wcj4kigDnezhLP48zSTuY9POL4vOfv71OBoEds4vQ5R0aJFefbZZy87PmDAAJcUJCKFTF66r3Lr4vV5HI6MT6dy2gXVqRPcfXfO1zC6pP7SHGcx9zKRZxjEOBbRma3U5eOvf6JBiyvU3iK7BiLiDk4HotmzZ2f7+UMu3eBHRDxCixYwfHjO2l3KFTOocuPCU57M1iHK6fo8zoSTTOp3AP2ZSGM20IWP2UM1mr5WmTElTe9ZrhdyFBGXc3phxtKlS2d4f+7cOU6fPo2fnx/Fixfn6NGjLi3QE2hhRvF6qakQGpr9OKKQEDh0KNOxNVSpYgYWZ7WacmYbl7pKblaqzu11srnP4wTzWMAcFv5zJwAdOpjtP8qUcX0pImK4dWHGY8eOZXidPHmS3bt307RpU+bNm5frokXEg/n6wrRp2beZNi3zoGH3DKr8Wp/nCvcZ7Ejikw/PMGUK+PnBF19AnTqwcaN7yhER5zgdiDJz7bXXMnr0aPr37++K04mIJ+rUCRYuzHyg8cKF2XdBecsMqivcp6NzJ556Cr79FqpXN2PJmzWDsWMhLc2ekkXEcMleZgAJCQk0b96c5ORkV5zOo6jLTOQieemCyq/uK7vl4D6Tk+Hxx2H+fPO+XTuYPdusZyQiruHM97fTgejzzz/P8N6yLA4ePMjkyZMJDw9n2bJlzlfs4RSIRMQdLAvefx+eeQbOnDEPlubN04LUIq7i1kDk45Oxl83hcFC2bFluu+02xo8fT3lXzxTxAApEIuJO330HDzwAu3eDjw+89prZI83HJYMaRLyXWwORN1IgEhF3O3kSnnoKPvzQvL/9dvPn0FB76xIpyNw6y+xilmWhPCUiknclS5oxRDNmQEAAxMSYWWirV9tdmYh3yFUgmj17NrVr1yYgIICAgABuuukmPrzw1xoREcm1hx+GLVvghhsgMRFat4ZXX3XtVm8icjmnA9GECRN48sknufPOO1mwYAELFizgjjvu4IknnuDNN990R40iIl6lVi3YtAl69zYDr4cNM11ortzdREQycnoMUdWqVRk2bNhlW3TMmjWLV199lT179ri0QE+gMUQiYpePPoInnoBTp6BcOfP+9tvtrkqkYHDrGKKDBw/SuHHjy443btyYg/rri4iIS/XoAfHxcNNNcPgwtG0LL78M58/bXZlI4eJ0IKpevToLFiy47PjHH3/Mtdde65KiRETkf66/3qxu/fjjpgttxAi47baMe9aKSN44vdv9sGHD6NKlC+vWraNJkyYAfPPNN6xatSrToCQiInkXEABTp0LLltCnj1kIu04dMzPtzjvtrk6k4HP6CVHnzp2Ji4vjqquuYsmSJSxZsoSrrrqKTZs2ce+997qjRhER+VeXLrB1K9StC0eOQPv28NxzcO6c3ZWJFGxamDEHNKhaRDxNSgoMGgSTJpn3kZFmX7Srr7a3LhFP4sz3d467zHK6aasCg4iI+/n7w8SJ0KKFmZ6/caPpQps5Ezp2tLk4kQIox0+IfHx8cDgcWX5uWRYOh4PUQrh6mJ4QiYgn27PHdKVt3mzeR0fDmDHg52drWSK2c8sTotjY2PQ/W5bFnXfeyfvvv0/FihVzX6mIiORZ1aqwfr3ZEHbCBHjrLfjmG/j4Y/OZiFxZrscQlSpViu3bt1OtWjVX1+Rx9IRIRAqKL76AXr3g2DEICoL//hc6dbK7KhF75NvmriIi4lk6dICEBGjcGJKSoHNnePppOHPG7spEPJsCkYhIIXP11bBmDTz/vHk/ebIJSL/8YmtZIh4tT4Eou0HWIiJin6JFYfRoWLYMrroKtm0zaxd9/LHdlYl4phwPqu50SSf0mTNneOKJJyhRokSG44sWLXJNZSIikmd33GG60B58ENatg65dITYW3nzTrH4tIkaOnxAFBQVlePXo0YMKFSpcdlxERDxLxYqwapXZFNbhgPfeg0aNYPduuysT8RxaqToHNMtMRAqLmBjo0QMOH4YSJcz+aD162F2ViHtolpmIiGTq9ttNF1rLlnDqFPTsaVa6Pn3a7spE7KVAJCLiZcqXN0+Khg0DHx+YMQNuuQV27bK7MhH7KBCJiHghX18YMsSMLQoLg++/N6FoxgzQQArxRgpEIiJerEUL2L4d2rSBf/4x3We9esHJk3ZXJpK/FIhERLxcuXJmvaIRI0wX2ocfQv368N13dlcmkn8UiEREBB8fePFFs8J1xYpmSn7DhjBtmrrQxDsoEImISLpmzcwstDvvNPufPf64WdQxOdnuykTcS4FIREQyuOoq+OILGDsWihSB+fOhXj2z/YdIYaVAJCIil/HxgUGDzHYfV19tNoZt1AimTFEXmhROCkQiIpKlyEjzZKhjRzh7Fvr1g/vvh+PH7a5MxLUUiEREJFtlysCSJWZD2KJFYeFCqFsXNm+2uzIR11EgEhGRK3I4IDoavvkGqlaFPXugSRN46y11oUnhoEAkIiI5dsstsHUrdO4M587BgAFwzz1w9KjdlYnkjQKRiIg4JTgYPvkEJk8GPz/4/HOIiICNG+2uTCT3FIhERMRpDgdERZkQdM01sHcvNG8O48ZBWprd1Yk4T4FIRERyrW5d04XWpQucPw/PPQcdOsDff9tdmYhzFIhERCRPAgNh3jx47z3w94f/+z+oUwe+/truykRyToFIRETyzOGAvn1h0ya4/nrYvx9atoSRI9WFJgWDApGIiLjMTTfBli3QowekpsJLL8Edd8Dhw3ZXJpI9BSIREXGpkiVh9mz4738hIABiYuDmmyE21u7KRLKmQCQiIi7ncMAjj5jVrGvVgsREaN0ahg0zT45EPI0CkYiIuM0NN5hxRY88YsYSvfoqtGkDBw/aXZlIRgpEIiLiViVKmO6z2bPNn1evNrPQYmLsrkzkfwpUIBo9ejQOh4Po6Oj0Y2fOnCEqKoqQkBBKlixJ586dOXToUIbf27t3L+3bt6d48eKUK1eOQYMGcf78+XyuXkTEu/XsaQZc165tBlm3bQsvv2zWLxKxW4EJRJs3b+a9997jpptuynB8wIABfPHFF3zyySesXbuWAwcO0KlTp/TPU1NTad++PWfPnmXDhg3MmjWLmTNnMmTIkPy+BRERr1ejBsTFmSn6lgUjRsBtt5lp+iJ2KhCB6OTJk3Tv3p3p06dTunTp9ONJSUl88MEHTJgwgdtuu4169eoxY8YMNmzYwLfffgvAV199xffff89HH31EnTp1aNeuHa+//jpTpkzh7NmzmV4vJSWF5OTkDC8REXGNgACziOO8eVCqlFnAsU4dWLbM7srEmxWIQBQVFUX79u1p3bp1huPx8fGcO3cuw/EaNWpw9dVXs/HfXQY3btxI7dq1CQ0NTW/Ttm1bkpOT2bVrV6bXGzVqFEFBQemv8PBwN9yViIh369oV4uPNxrB//w133gnPPw/nztldmXgjjw9E8+fPZ+vWrYwaNeqyzxITE/Hz8yM4ODjD8dDQUBITE9PbXByGLnx+4bPMDB48mKSkpPTXvn37XHAnIiJyqWuvhQ0bzEaxAGPHQosWZrNYkfzk0YFo37599O/fnzlz5lCsWLF8u66/vz+BgYEZXiIi4h7FisHkyfDppxAUZAJSnTrw+ed2VybexKMDUXx8PIcPH6Zu3boUKVKEIkWKsHbtWiZOnEiRIkUIDQ3l7NmzHD9+PMPvHTp0iLCwMADCwsIum3V24f2FNiIiYr/OnWHrVrjlFjh2DO6+GwYOhCyGe4q4lEcHolatWrFjxw4SEhLSX/Xr16d79+7pfy5atCirVq1K/53du3ezd+9eIiMjAYiMjGTHjh0cvmgjnZiYGAIDA6lVq1a+35OIiGStWjVYvx4urK7y5pvQrBns2WNrWeIFithdQHZKlSrFjTfemOFYiRIlCAkJST/+6KOPMnDgQMqUKUNgYCBPP/00kZGRNGrUCIA2bdpQq1YtevbsydixY0lMTOTll18mKioKf3//fL8nERHJnp+fCUItW8LDD5uVriMizOKOF62qIuJSHv2EKCfefPNN7rrrLjp37kzz5s0JCwtj0aJF6Z/7+vqydOlSfH19iYyMpEePHjz00EO89tprNlYtIiJX0rEjbNsGjRpBUpLpUnv6aUhJsbsyKYwclmVZdhfh6ZKTkwkKCiIpKUkDrEVE8tm5c2ZF67Fjzfu6deHjj6F6dXvrEs/nzPd3gX9CJCIihVvRojBmDHz5JYSEmIHXF0KRiKsoEImISIFw552QkABNm8KJE2ZhxyeegH/+sbsyKQwUiEREpMCoVAliY+Gll8DhMFuANGoEu3fbXZkUdApEIiJSoBQpAsOHw4oVULYsfPcd1KsHH31kd2VSkCkQiYhIgXT77bB9u5mef+oU9OwJjz4Kp0/bXZkURApEIiJSYJUvDzEx8Oqrpgvtv/+FBg3g++/trkwKGgUiEREp0Hx9YehQWLUKwsJg1y6oXx9mzrS7MilIFIhERKRQaNnSzEK7/XYz8+yRR6BXLzh50u7KpCBQIBIRkUIjNBSWL4cRI8DHB2bPNpvFfved3ZWJp1MgEhGRQsXHB158EdasgYoV4ccfoWFDmDYNtDeDZEWBSERECqVmzUwXWrt2cOYMPP44PPggJCfbXZl4IgUiEREptK66CpYuNfug+frC/PlmzaJt2+yuTDyNApGIiBRqPj4waBB8/TWEh8Mvv5jVrd95R11o8j8KRCIi4hUiI00XWseOcPYsREXBAw9AUpLdlYknUCASERGvUaYMLFkCEyZA0aLw6acQEQGbN9tdmdhNgUhERLyKwwEDBsD69VClCuzZA02awNtvqwvNmykQiYiIV2rQwAyu7tQJzp2D6Gi49144etTuysQOCkQiIuK1goNNt9mkSeDnB599ZrrQvv3W7sokvykQiYiIV3M4oF8/2LgRrrkG9u41axiNGwdpaXZXJ/lFgUhERASoWxe2boUuXeD8eXjuOTMj7e+/7a5M8oMCkYiIyL8CA2HePJg6Ffz94csvoU4dMwBbCjcFIhERkYs4HGabj7g4uO462L8fWrSAUaPUhVaYKRCJiIhk4uabIT4eevSA1FSzYWy7dnD4sN2ViTsoEImIiGShZEmYPRs++AACAuCrr0wX2po1dlcmrqZAJCIikg2HA3r3hk2boGZNOHgQWrWC114zT46kcFAgEhERyYEbbzRbfDzyiBlLNHQotGkDiYl2VyauoEAkIiKSQyVKwH//a7rRiheH1avNWKOVK+2uTPJKgUhERMRJPXuaAde1a5tB1m3awCuvmPWLpGBSIBIREcmFGjXM1Py+fc2msMOHm7FF+/fbXZnkhgKRiIhILgUEwHvvmcUcS5aEdevMLLTly+2uTJylQCQiIpJHXbuabT8iIsxWH+3awQsvwLlzdlcmOaVAJCIi4gLXXgsbNkBUlHk/ZoxZ4XrvXlvLkhxSIBIREXGRYsVg8mT49FMICjIBKSICvvjC7srkShSIREREXKxzZ9OFdsstcPQodOwI//kPnD1rd2WSFQUiERERN6hWDdavh+ho837CBGjWDPbssbUsyYICkYiIiJv4+cGbb8KSJRAcbLb/iIiAxYvtrkwupUAkIiLiZnffDQkJ0KgRJCVBp07wzDOQkmJ3ZXKBApGIiEg+qFzZrFM0aJB5P2kSNG4Mv/xib11iKBCJiIjkk6JFYexYWLoUQkLMwOu6dWHBArsrEwUiERGRfNa+velCa9oUTpyALl3gySfhn3/srsx7KRCJiIjYoFIliI2FF18EhwOmTjVjjHbvtrsy76RAJCIiYpMiRWDECLP3Wdmy8N13UK8ezJljd2XeR4FIRETEZm3awPbtZquPU6egRw947DE4fdruyryHApGIiIgHKF8eVq6EoUNNF9oHH0CDBvD993ZX5h0UiERERDyEry+8+qoJRmFhsGuX2f5j5ky7Kyv8FIhEREQ8zG23mVlorVubbrNHHoFeveDkSbsrK7wUiERERDxQaCisWAHDh4OPD8yebZ4W7dhhd2WFkwKRiIiIh/LxgZdeMtPzK1SAH38044qmTwfLsru6wkWBSERExMM1b2660O64A86cgb59oXt3s6ijuIYCkYiISAFQtix8+SWMGWMGX8+bZ7b92LbN7soKBwUiERGRAsLHB557zmwSGx5uNoaNjIR33lEXWl4pEImIiBQwjRubJ0MdOkBKCkRFwQMPQFKS3ZUVXApEIiIiBVBICHz2GYwfb7YA+fRT04W2ZYvdlRVMCkQiIiIFlMMBAwfC+vVQuTL89pt5evT22+pCc5YCkYiISAHXsKHpQrv3Xjh3DqKjzZ+PHbO7soJDgUhERKQQKF0aFi6EiRPBz890p0VEwLff2l1ZwaBAJCIiUkg4HPD007BhA1xzDfzxBzRrBm+8AWlpdlfn2RSIRERECpl69SA+3sw8O38eBg2Cjh3hyBG7K/NcCkQiIiKFUFAQzJ8P774L/v5mUcc6dcwAbLmcApGIiEgh5XDAE09AXBxcdx38+Se0aAGjRqkL7VIeHYhGjRrFLbfcQqlSpShXrhz33HMPu3fvztDmzJkzREVFERISQsmSJencuTOHDh3K0Gbv3r20b9+e4sWLU65cOQYNGsT58+fz81ZERERsc/PNZn2i7t0hNRVefBHuvBMOH7a7Ms/h0YFo7dq1REVF8e233xITE8O5c+do06YNp06dSm8zYMAAvvjiCz755BPWrl3LgQMH6NSpU/rnqamptG/fnrNnz7JhwwZmzZrFzJkzGTJkiB23JCIiYotSpeDDD+GDDyAgAFasMF1oa9faXZlncFhWwVm66a+//qJcuXKsXbuW5s2bk5SURNmyZZk7dy733XcfAD/++CM1a9Zk48aNNGrUiGXLlnHXXXdx4MABQkNDAZg6dSrPP/88f/31F35+fle8bnJyMkFBQSQlJREYGOjWexQREXG3nTvNgOsffjD7ow0dCi+9ZDaNLUyc+f726CdEl0r6d5OWMmXKABAfH8+5c+do3bp1epsaNWpw9dVXs3HjRgA2btxI7dq108MQQNu2bUlOTmbXrl2ZXiclJYXk5OQMLxERkcLixhth82Z45BEzlmjoUGjbFhIT7a7MPgUmEKWlpREdHU2TJk248cYbAUhMTMTPz4/g4OAMbUNDQ0n8959qYmJihjB04fMLn2Vm1KhRBAUFpb/Cw8NdfDciIiL2KlEC/vtfmD0biheHVatMF9rKlXZXZo8CE4iioqLYuXMn8+fPd/u1Bg8eTFJSUvpr3759br+miIiIHXr2NGsW1a4Nhw5Bmzbwyitm/SJvUiACUb9+/Vi6dCmxsbFUqlQp/XhYWBhnz57l+PHjGdofOnSIsLCw9DaXzjq78P5Cm0v5+/sTGBiY4SUiIlJY1ahhpub37Ws2hR0+HFq1gv377a4s/3h0ILIsi379+rF48WJWr15N1apVM3xer149ihYtyqpVq9KP7d69m7179xIZGQlAZGQkO3bs4PBFcwtjYmIIDAykVq1a+XMjIiIiHi4gAN57D+bNg5IlYd0604W2fLndleUPj55l9tRTTzF37lw+++wzrr/++vTjQUFBBAQEAPDkk0/yf//3f8ycOZPAwECefvppADZs2ACYafd16tShQoUKjB07lsTERHr27Mljjz3GyJEjc1SHZpmJiIg3+flnMwstIcG8f+EFeP11KFLE1rKc5sz3t0cHIofDkenxGTNm8PDDDwNmYcb//Oc/zJs3j5SUFNq2bcs777yToTvsjz/+4Mknn2TNmjWUKFGCXr16MXr0aIrk8J+sApGIiHibM2fg2WdhyhTzvkkT8/SoIM0zKjSByFMoEImIiLf69FN49FFIToYyZWDWLLjrLruryplCuw6RiIiI5K/77oNt26B+fTh6FDp0ME+Ozp61uzLXUiASERGRbFWrBuvXQ3S0eT9+PDRvDr//bmdVrqVAJCIiIlfk7w9vvglLlkBwsJmmHxFh3hcGCkQiIiKSY3ffbWafNWoEx4/DvfdC//6QkmJ3ZXmjQCQiIiJOqVzZrFM0aJB5P3GimYX266/21pUXCkQiIiLitKJFYexYWLoUQkLM9h9168Inn9hdWe4oEImIiEiutW9vutCaNjVT8x94AJ56yqxjVJAoEImIiEieVKoEsbEweLB5/+67EBlpVrwuKBSIREREJM+KFIGRI83eZ2XLmqdGdeua1a0LAgUiERERcZm2bU0YatECTp6EBx+EPn3g9Gm7K8ueApGIiIi4VIUKsHIlDBkCDge8/z40bAg//GB3ZVlTIBIRERGX8/WFYcNMMAoNhZ07zfYfs2bZXVnmFIhERETEbW67DbZvh9atTbfZww+b16lTdleWkQKRiIiIuFVoqBls/frr4ONjnhLdcot5auQpFIhERETE7Xx94eWXYfVqM8bohx9MKHr/fbAsu6tTIBIREZF8dOutZhbaHXeYxRv79IEePeDECXvrUiASERGRfFW2LHz5JYwebZ4czZ0L9erBwYP21aRAJCIiIvnOxweefx7WrjUrXVerZsYa2aWIfZcWERERb9ekielCsywTkuyiQCQiIiK2CgmxuwJ1mYmIiIgoEImIiIgoEImIiIjXUyASERERr6dAJCIiIl5PgUhERES8ngKRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEImIiIjXUyASERERr6dAJCIiIl5Pu93ngGVZACQnJ9tciYiIiOTUhe/tC9/j2VEgyoETJ04AEB4ebnMlIiIi4qwTJ04QFBSUbRuHlZPY5OXS0tI4cOAApUqVwuFwuPTcycnJhIeHs2/fPgIDA116bk+i+yxcvOU+wXvuVfdZuOg+DcuyOHHiBBUqVMDHJ/tRQnpClAM+Pj5UqlTJrdcIDAws1P+nvUD3Wbh4y32C99yr7rNw0X1yxSdDF2hQtYiIiHg9BSIRERHxegpENvP392fo0KH4+/vbXYpb6T4LF2+5T/Cee9V9Fi66T+dpULWIiIh4PT0hEhEREa+nQCQiIiJeT4FIREREvJ4CkYiIiHg9BSIbTZkyhSpVqlCsWDEaNmzIpk2b7C7J5datW0eHDh2oUKECDoeDJUuW2F2SW4waNYpbbrmFUqVKUa5cOe655x52795td1ku9+6773LTTTelL4IWGRnJsmXL7C7L7UaPHo3D4SA6OtruUlzu1VdfxeFwZHjVqFHD7rLcYv/+/fTo0YOQkBACAgKoXbs2W7Zssbssl6pSpcpl/zwdDgdRUVF2l+ZSqampvPLKK1StWpWAgACuueYaXn/99RztWZYVBSKbfPzxxwwcOJChQ4eydetWbr75Ztq2bcvhw4ftLs2lTp06xc0338yUKVPsLsWt1q5dS1RUFN9++y0xMTGcO3eONm3acOrUKbtLc6lKlSoxevRo4uPj2bJlC7fddht33303u3btsrs0t9m8eTPvvfceN910k92luM0NN9zAwYMH01/r16+3uySXO3bsGE2aNKFo0aIsW7aM77//nvHjx1O6dGm7S3OpzZs3Z/hnGRMTA8D9999vc2WuNWbMGN59910mT57MDz/8wJgxYxg7diyTJk3K/UktsUWDBg2sqKio9PepqalWhQoVrFGjRtlYlXsB1uLFi+0uI18cPnzYAqy1a9faXYrblS5d2nr//fftLsMtTpw4YV177bVWTEyMdeutt1r9+/e3uySXGzp0qHXzzTfbXYbbPf/881bTpk3tLiPf9e/f37rmmmustLQ0u0txqfbt21u9e/fOcKxTp05W9+7dc31OPSGywdmzZ4mPj6d169bpx3x8fGjdujUbN260sTJxlaSkJADKlCljcyXuk5qayvz58zl16hSRkZF2l+MWUVFRtG/fPsO/q4XRzz//TIUKFahWrRrdu3dn7969dpfkcp9//jn169fn/vvvp1y5ckRERDB9+nS7y3Krs2fP8tFHH9G7d2+Xb0xut8aNG7Nq1Sp++uknALZv38769etp165drs+pzV1t8Pfff5OamkpoaGiG46Ghofz44482VSWukpaWRnR0NE2aNOHGG2+0uxyX27FjB5GRkZw5c4aSJUuyePFiatWqZXdZLjd//ny2bt3K5s2b7S7FrRo2bMjMmTO5/vrrOXjwIMOGDaNZs2bs3LmTUqVK2V2ey/z222+8++67DBw4kBdffJHNmzfzzDPP4OfnR69evewuzy2WLFnC8ePHefjhh+0uxeVeeOEFkpOTqVGjBr6+vqSmpjJixAi6d++e63MqEIm4WFRUFDt37iyU4zAArr/+ehISEkhKSuLTTz+lV69erF27tlCFon379tG/f39iYmIoVqyY3eW41cV/o77pppto2LAhlStXZsGCBTz66KM2VuZaaWlp1K9fn5EjRwIQERHBzp07mTp1aqENRB988AHt2rWjQoUKdpficgsWLGDOnDnMnTuXG264gYSEBKKjo6lQoUKu/3kqENngqquuwtfXl0OHDmU4fujQIcLCwmyqSlyhX79+LF26lHXr1lGpUiW7y3ELPz8/qlevDkC9evXYvHkzb7/9Nu+9957NlblOfHw8hw8fpm7duunHUlNTWbduHZMnTyYlJQVfX18bK3Sf4OBgrrvuOn755Re7S3Gp8uXLXxbaa9asycKFC22qyL3++OMPVq5cyaJFi+wuxS0GDRrECy+8QNeuXQGoXbs2f/zxB6NGjcp1INIYIhv4+flRr149Vq1alX4sLS2NVatWFdqxGIWdZVn069ePxYsXs3r1aqpWrWp3SfkmLS2NlJQUu8twqVatWrFjxw4SEhLSX/Xr16d79+4kJCQU2jAEcPLkSX799VfKly9vdyku1aRJk8uWwvjpp5+oXLmyTRW514wZMyhXrhzt27e3uxS3OH36ND4+GSOMr68vaWlpuT6nnhDZZODAgfTq1Yv69evToEED3nrrLU6dOsUjjzxid2kudfLkyQx/09yzZw8JCQmUKVOGq6++2sbKXCsqKoq5c+fy2WefUapUKRITEwEICgoiICDA5upcZ/DgwbRr146rr76aEydOMHfuXNasWcOKFSvsLs2lSpUqddn4rxIlShASElLoxoU9++yzdOjQgcqVK3PgwAGGDh2Kr68v3bp1s7s0lxowYACNGzdm5MiRPPDAA2zatIlp06Yxbdo0u0tzubS0NGbMmEGvXr0oUqRwfs136NCBESNGcPXVV3PDDTewbds2JkyYQO/evXN/0jzOfJM8mDRpknX11Vdbfn5+VoMGDaxvv/3W7pJcLjY21gIue/Xq1cvu0lwqs3sErBkzZthdmkv17t3bqly5suXn52eVLVvWatWqlfXVV1/ZXVa+KKzT7rt06WKVL1/e8vPzsypWrGh16dLF+uWXX+wuyy2++OIL68Ybb7T8/f2tGjVqWNOmTbO7JLdYsWKFBVi7d++2uxS3SU5Otvr3729dffXVVrFixaxq1apZL730kpWSkpLrczosKw/LOoqIiIgUAhpDJCIiIl5PgUhERES8ngKRiIiIeD0FIhEREfF6CkQiIiLi9RSIRERExOspEImIiIjXUyASERERr6dAJCKF2po1a3A4HBw/ftzuUkTEgykQiYitEhMTefrpp6lWrRr+/v6Eh4fToUOHDJsf28HhcKS/AgMDueWWW/jss8+cOsfvv/+Ow+EgISHBPUWKiMsoEImIbX7//Xfq1avH6tWrGTduHDt27GD58uW0bNmSqKgou8tjxowZHDx4kC1bttCkSRPuu+8+duzYYXdZIuIGCkQiYpunnnoKh8PBpk2b6Ny5M9dddx033HADAwcO5Ntvv01vN2HCBGrXrk2JEiUIDw/nqaee4uTJk+mf//HHH3To0IHSpUtTokQJbrjhBv7v//4vw7Xi4+OpX78+xYsXp3HjxuzevfuK9QUHBxMWFsZ1113H66+/zvnz54mNjU3/fPny5TRt2pTg4GBCQkK46667+PXXX9M/r1q1KgARERE4HA5atGiR/tn7779PzZo1KVasGDVq1OCdd95x+n8/EXEdBSIRscXRo0dZvnw5UVFRlChR4rLPg4OD0//s4+PDxIkT2bVrF7NmzWL16tU899xz6Z9HRUWRkpLCunXr2LFjB2PGjKFkyZIZzvfSSy8xfvx4tmzZQpEiRejdu3eOaz1//jwffPABAH5+funHT506xcCBA9myZQurVq3Cx8eHe++9l7S0NAA2bdoEwMqVKzl48CCLFi0CYM6cOQwZMoQRI0bwww8/MHLkSF555RVmzZqV45pExMUsEREbxMXFWYC1aNEip3/3k08+sUJCQtLf165d23r11VczbRsbG2sB1sqVK9OPffnllxZg/fPPP1leA7CKFStmlShRwvLx8bEAq0qVKtaRI0ey/J2//vrLAqwdO3ZYlmVZe/bssQBr27ZtGdpdc8011ty5czMce/31163IyMgszy0i7qUnRCJiC8uyctx25cqVtGrViooVK1KqVCl69uzJkSNHOH36NADPPPMMw4cPp0mTJgwdOpTvvvvusnPcdNNN6X8uX748AIcPH872um+++SYJCQksW7aMWrVq8f7771OmTJn0z3/++We6detGtWrVCAwMpEqVKgDs3bs3y3OeOnWKX3/9lUcffZSSJUumv4YPH56hu01E8pcCkYjY4tprr8XhcPDjjz9m2+7333/nrrvu4qabbmLhwoXEx8czZcoUAM6ePQvAY489xm+//UbPnj3ZsWMH9evXZ9KkSRnOU7Ro0fQ/OxwOgPSurayEhYVRvXp12rRpw4wZM+jSpUuGENWhQweOHj3K9OnTiYuLIy4uLkNdmbkw9mn69OkkJCSkv3bu3Jlh3JSI5C8FIhGxRZkyZWjbti1Tpkzh1KlTl31+Yd2g+Ph40tLSGD9+PI0aNeK6667jwIEDl7UPDw/niSeeYNGiRfznP/9h+vTpLq23QYMG1KtXjxEjRgBw5MgRdu/ezcsvv0yrVq2oWbMmx44dy/A7F8Ybpaamph8LDQ2lQoUK/Pbbb1SvXj3D68IgbBHJfwpEImKbKVOmkJqaSoMGDVi4cCE///wzP/zwAxMnTiQyMhKA6tWrc+7cOSZNmsRvv/3Ghx9+yNSpUzOcJzo6mhUrVrBnzx62bt1KbGwsNWvWdHm90dHRvPfee+zfv5/SpUsTEhLCtGnT+OWXX1i9ejUDBw7M0L5cuXIEBASwfPlyDh06RFJSEgDDhg1j1KhRTJw4kZ9++okdO3YwY8YMJkyY4PKaRSSH7B7EJCLe7cCBA1ZUVJRVuXJly8/Pz6pYsaLVsWNHKzY2Nr3NhAkTrPLly1sBAQFW27ZtrdmzZ1uAdezYMcuyLKtfv37WNddcY/n7+1tly5a1evbsaf3999+WZf1vUPWFtpZlWdu2bbMAa8+ePVnWBViLFy/OcCwtLc2qUaOG9eSTT1qWZVkxMTFWzZo1LX9/f+umm26y1qxZc9nvTZ8+3QoPD7d8fHysW2+9Nf34nDlzrDp16lh+fn5W6dKlrebNm+dqgLmIuIbDspwY2SgiIiJSCKnLTERERLyeApGIiIh4PQUiERER8XoKRCIiIuL1FIhERETE6ykQiYiIiNdTIBIRERGvp0AkIiIiXk+BSERERLyeApGIiIh4PQUiERER8Xr/D9pZay0zUGKXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "keep = cash_rate < 6\n",
    "cash_rate = cash_rate[keep]\n",
    "house_prices = house_prices[keep]\n",
    "\n",
    "plt.plot(cash_rate, house_prices, 'ro')\n",
    "plt.xlabel(\"Cash Rate\")\n",
    "plt.ylabel(\"House Prices\")\n",
    "\n",
    "ones = np.ones((len(cash_rate)))\n",
    "X = np.vstack([ones, cash_rate]).T\n",
    "Y = house_prices.reshape((len(house_prices), 1))\n",
    "beta = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(Y))\n",
    "x_model = np.arange(0, 8, 0.25)\n",
    "y_model = x_model * beta[1] + beta[0]  # Use our learned model to fit the line of best fit.\n",
    "\n",
    "plt.plot(x_model, y_model, 'b-')\n",
    "predicted_prices = cash_rate * beta[1] + beta[0]\n",
    "squared_error = np.sum((predicted_prices - house_prices) ** 2)\n",
    "print(\"The squared error is {:.2f}\".format(squared_error))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1072.17295516]\n",
      " [-123.36619231]]\n"
     ]
    }
   ],
   "source": [
    "# Question 2\n",
    "\n",
    "# Store parameter versions here\n",
    "models = []\n",
    "\n",
    "n = len(cash_rate)\n",
    "\n",
    "for i in range(100):\n",
    "    \n",
    "    # Choose random indexes\n",
    "    # replace=True is better statistically, we will come back to it later on...\n",
    "    random_indexes = np.random.choice(np.arange(n), size=int(n/2), replace=True)\n",
    "    cash_rate_sample = cash_rate[random_indexes]\n",
    "    house_prices_sample = house_prices[random_indexes]\n",
    "    \n",
    "    ones_sample = np.ones((len(cash_rate_sample)))\n",
    "    X_sample = np.vstack([ones_sample, cash_rate_sample]).T\n",
    "\n",
    "    Y_sample = house_prices_sample.reshape((len(house_prices_sample), 1))\n",
    "\n",
    "    beta_sample = np.linalg.inv(X_sample.T.dot(X_sample)).dot(X_sample.T.dot(Y_sample))\n",
    "    models.append(beta_sample)\n",
    "# Average all parameters\n",
    "final_model = np.mean(models, axis=0)\n",
    "print(final_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/ols_basic_ensemble.py`*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
